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SUMMARY 


Our  analysis  of  the  free  sheared  atmosphere  is  based  on  the  hierarchy  of  equations  satisfied 
by  the  modal  coefficients  (strengths).  This  hierarchy  is  obtained  by  substituting  a  series 
expansion  of  the  fluid  variables  into  the  governing  equations.  Our  analysis  is  particularly 
focused  on  the  stability  of  the  Taylor-Dyson  atmosphere,  which  is  designed  to  be  a  model  of 
the  free  sheared  atmosphere.  The  construction  of  the  hierarchy  is  however  independent  of 
the  model  that  one  chooses  albeit  the  more  complicated  the  model,  the  more  complicated  the 
hierarchy.  The  classical  approach  considers  Fourier  scries  and  the  prototype  is  the  successful 
Salzman- Lorenz  hierarchy  for  the  Benard  problem,  i.e.  a  thin  layer  of  fluid  heated  from 
below.  This  problem  is  designed  to  be  a  model  of  the  lower  atmosphere.  The  application 
to  the  free  sheared  atmosphere  modeled  by  the  famous  Taylor-Dyson  equations  of  Fourier 
techniques  has  yielded  the  puzzling  result  that  no  unstable  modes  appear  in  linear  order.  The 
origin  of  the  Richardson  number  criterion  thus  remains  obscure.  With  the  advent  of  wavelet 
analysis,  we  have  deemed  it  useful  (actually  imperative)  to  perform  the  analysis,  not  on  the 
basis  of  the  Fourier  modes  which  have  infinite  energy  and  essentially  unbounded  support, 
but  rather  on  the  basR  ot  wavelets  which  are  designed  to  have  finite  energy  and  essentially 
compact  support.  It  is  intuitively  clear  that  such  refined  modal  analysis  should  allow  for 
a  much  more  satisfactory  physical  interpretation  of  results.  This  point  of  view  is  further 
strengthened  by  the  current  trend  in  experimental  work  that  emphasizes  disturbances  with 
finite  spatial  and  temporal  extents.  At  Boston  University,  we  have  developed  a  full  form  of 
wavelet  expansion  which  has  the  advantage  over  more  conventional  procedures  of  possessing 
a  complete  convolution  algebra. 

The  possibility  of  operating  with  a  convolution  algebra  has  so  far  been  a  prerogative  of  the 
Fourier  analysis.  Its  major  use  in  the  present  context  is  the  calculation  of  the  effect  of  forcing 
on  a  particular  hierarchical  level.  Convolution  algebra  permits  an  immediate  expression  for 
the  result  of  a  forcing  term  on  a  variable  that  satisfies  a  given  differential  equation.  The 


IX 


mathematical  >ols  needed  to  obtain  a  convolution  algebra  for  modes  with  finite  energy 
are  rather  extensive  and  not  always  elementary.  To  help  our  thinking,  and  also  to  help 
introduce  interested  scientist  to  a  rather  difficult  area,  we  put  together  the  most  essential 
parts  of  the  analysis  and  we  present  them  in  this  report.  It  is  worth  mentioning  at  this  point 
that  very  fruitful  applications  have  been  made  of  our  wavelet  analysis  in  biomedical  areas. 
In  particular  the  analysis  of  M  waves  in  muscle  response  has  been  analysed  very  successfully 
with  the  wavelets  presented  here.  While  our  main  interest  is  the  development  of  the  chaos 
dynamics  of  the  free  sheared  atmosphere,  the  Taylor-Dyson  model  in  particular,  we  deem 
the  mathematical  background  essential  and  also  of  general  interest. 
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1  Signals  and  Systems  and  the  Inverse  Problem 


1.1  Introduction 

A  majority  of  systems,  natural  and  man-made,  are  linear  and  time-invariant  ( LT I ) [  1 J . 
By  linear,  it  is  meant  that  if  ut  and  u2  are  two  inputs  to  a  physical  system,  with 
respective  outputs  yx  and  y2,  then  the  output  to  the  input  aux  -f  0u2  is  aj/!  +  0y2, 
where  n  and  0  arc  two  real  numbers.  By  time-invariant,  it  is  meant  that  a  time- 
shift  in  the  input  signal  causes  an  equal  time-shift  in  the  output  signal.  In  other 
words,  if  the  input  u(t)  produces  the  output  y{t),  then  the  input  u(t  —  r)  produces 
the  output  y(t  —  r).  These  properties  are  depicted  in  figure  1.  It  can  be  proven  that 
the  output  of  any  linear  time-invariant  system  can  be  modelled  by  the  convolution  of 
the  input  signal  with  the  impulse  response  of  that  system.  The  impulse  response  (or 
point-spread  function  PSF)  of  a  system  is  the  output  when  the  input  to  the  system 
is  a  Dirac  delta  function  (e.g.  for  electrical  systems,  a  very  high,  very  short-duration 
electrical  pulse).  The  proof  is  as  follows.  Let  C  denote  a  discrete  LTI  system.  Let 
a  discrete  Dirac  delta  <5[n]  be  the  input  to  the  system.  By  definition,  the  discrete 
impulse  response  h0[n}  is  the  output  of: 

<5[n]  — ►  [Z]  — ♦  h0[n\  (1) 

Because  of  the  property  of  shift-invariance, 

S[n  -  k]  — >  [Zj — ♦  hk[n]  =  h0[n  -  k]  (2) 

Furthermore,  because  of  the  property  of  linearity, 

x[£]<5[n  -  fc]  — ♦  [Z]  — ♦  XZMM"]  v3) 

k  k 

It  can  be  proven  that  any  discrete  signal  x[n]  can  be  written  in  the  form  x[u]  = 
Efc  •r[k]<5[n  -  k}.  Th  is  property  is  commonly  referred  to  as  the  discrete  sampling 
propertxj[  l],  or,  more  technically,  as  Dirac's  identity.  Combining  Eqs.(2)  and  (3),  we 

obtain  _ 

x[n]  =  ]Z  x[A:]c!)[n  -  A;]  - ►  C_  - ♦  y[n]  =  ^  j[A-]/?0[??  —  k]  (-1) 

k  k 

In  the  case  of  continuous-time  LTI  systems,  Eq.('l)  is  written  as: 

/  +  OO  -  r  +  OO 

x(t  )S(t  —  t)  dr  - >  [Zj  - »  y(t)  =  /  x(r)/i0(/ -r )  dr  =  x(t)  *  h0(  0 

-  OO  J  —  OO 

(5) 

where  /i0(/)  is  the  impulse  response  of  C. 

In  many  instances,  the  physicist,  the  engineer,  or  the  mathematician  is  asked  to 
determine  the  initial  condition  (or  input)  to  a  system,  when  the  final  condition  (or 


output)  is  known.  This  problem  is  known  as  the  inverse  problem.  In  mathematical 
terminology,  knowledge  of  the  system  implies  knowledge  of  its  impulse  response;  the 
solution  of  the  inverse  problem  implies  knowledge  of  the  output,  and  convolution 
inverse  to  its  impulse  response.  This  last  statement  is  proven  below.  If  u  represents 
the  input  to  the  system,  y  the  output,  and  h  its  impulse  response,  then  by  definition: 

/+oo 

h(x')u{x  —  x')dx'  (£) 

-OO 

If  the  impulse  response  h  has  a  convolution  inverse  7nu[A],  this  implies  by  definition 
that  the  convolution  of  h  with  Inv[h]  yields  the  Dirac  delta  function: 

h(x)  *  7nt>[/i](:r)  =  6(x)  (7) 

That  is  because  the  Dirac  delta  is  the  unit  element  for  the  operation  of  convolution, 
much  in  the  same  way  the  number  0  is  the  unit  element  for  integer  addition: 

f(x)*6(x)  =  f(x)  for  any  function  f(x)  .  , 

n  +  0  =  n  for  any  integer  n 

If  y  an  2  h  are  known,  then  the  initial  condition  u  can  be  determined  in  the  following 
way:  convolve  Eq.(6)  with  /nu[A],  and  use  Eq.(7)  to  obtain: 

7nu[/i]  *  y  —  7nu[A]  *h*u  =  8*u  =  u  (9) 

With  the  input  u  thus  determined,  the  inverse  problem  for  a  linear  time-invariant 
system  can  always  be  solved,  provided  the  convolution  inverse  lnv[h]  to  the  system’s 
impulse  response  can  be  evaluated.  Unfortunately,  in  a  number  of  cases,  the  con¬ 
volution  inverse  cannot  be  determined  within  the  space  of  smooth  point-functions, 
much  in  the  same  way  natural  numbers  do  not  have  natural  number  inverses  to  the 
operation  of  addition:  one  either  has  to  construct  a  new  operation  (ic.  subtraction) 
or  extend  the  space  of  natural  numbers  to  the  space  or  positive  and  negative  integers: 


n  —  n  =  0,  or  (10) 

n  +  (-n)  =  0  (11) 

By  analogy,  the  construction  of  a  consistent  convolution  inverse  to  any  impulse  re¬ 
sponse  will  require: 


•  either  a  new  defining  operation,  replacing  convolution  as  the  operation  of  choice 
for  modelling  LTI  systems, 

•  or  the  extension  of  the  space  of  smooth  point-functions  to  a  larger  space,  which 
will  include  more  singular  functions  than  just  smooth  or  piecewise  smooth  point- 
functions. 
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This  picture  will  become  clearer  in  the  next  section,  when  Fourier  transforms  are 
introduced,  and  where  it  will  be  shown  that  there  exists  no  smooth  point-function  such 
that  convolving  with  any  Gaussian  e~x2^2  /  y/rtX  yields  a  Dirac  delta;  in  other  words, 
that  there  exists  no  smooth  or  piecewise  smooth  point-function  that  is  convolution 
inverse  to  a  Gaussian  (of  infinite  support,  i.e.  defined  in  ]  —  oo,  -foo[). 


1.2  Fourier  Transform  Solutions  to  the  Inverse  Problem 


Fourier  transforms  are  frequently  utilized  to  solve  inverse  problems  for  linear  time- 
invariant  systems.  That  is  because  the  operation  of  convolution  gets  mapped  into  a 
simple  algebraic  multiplication  under  the  Fourier  transform.  If  u(x),  y(x),  and  h(x ) 


are  respectively  the  input,  output,  and  impulse  response  of  a 
system,  and  if  y(k)  denotes  the  Fourier  transform  of  y(x): 

linear  time- invariant 

y{k)  =  f  y{x)e~'kz  dx 

J — OO 

(12) 

1  r+oo 

y(*)  =  ~  J  ^  y(k)e+tkxdk 

(13) 

then 

y(x )  =  u(x)  *  h(x) 

(H) 

y(k)  =  u(k)-h(k) 

(15) 

and  we  can  thus  obtain  an  expression  for  the  input,  given  the  output  and  the  impulse 
response: 


u(k) 

u(x) 


y(k ) 

h(k) 

±  r  ik 

2k  J-„  Mk) 


(16) 

(17) 


However,  eventual  zeros  of  h(k),  at  specific  points  in  wave-number  space,  or  at  the 
limits  ±oo  of  the  Fourier  integral,  will  ill-condition  the  inverse- Fourier  integral  in 
Eq.(17).  This  result  can  be  interpreted  in  the  following  way.  In  Fourier  space,  the 
blurred  waveform  is  the  product  of  the  optical  transfer  function  (OTF,  Fourier  trans¬ 
form  of  the  point-spread  function)  with  the  original  waveform.  As  the  support  of 
the  point-spread  function  increases,  the  support  of  the  optical  transfer  function  de¬ 
creases  (a  property  of  Fourier  transforms),  thus  cutting  ofT  an  increasing  amount  of 
high  frequencies  of  the  original  waveform,  which  unavoidably  “smooths”  or  “blurs” 
it.  This  process  is  depicted  in  figure  2.  The  inverse  Fourier  integral  in  Fq.(17)  has 
to  restore  the  high  frequencies  lost  in  the  blurring  process;  the  higher  the  frequen¬ 
cies,  the  greater  the  magnitude  of  the  restoration.  This  process  eventually  leads  to  a 
diverging  inverse- Fourier  integral. 


3 


We  now  consider  the  case  of  a  Gaussian  filter  G(x,t): 


G(x,t) 


_  il 

e 

•s/Att t 


(18) 


and  attempt  to  evaluate  an  explicit  representation  of  the  convolution  inverse  of  G(x,  t ) 
with  the  help  of  Fourier  transforms.  We  suppress  the  spatial  variable  x  if  confusion 
does  not  arise,  and  denote  G(x,  t)  by  G(t).  Thus,  using  standard  definitions[2][3],  we 
introduce 


F(a,  x,  2t) 


i  dk 

2x7-0 

£ 

^(erff.Vi  -  +  +  crf(is/i  a  -  ^ 


(19) 

(20) 


A  graph  of  F(a,x,l)  appears  in  figure  3,  and  the  time  history  of  F(a,x,t)  appears 
in  figure  4.  In  constructing  the  figures,  the  Nyquist  criterion  was  applied  to  ensure 
that  the  factor  did  not  generate  frequency  aliasing[l]. 

Proceeding  formally,  Inv[G(t)\  =  lima_00  F(a,x,2t),  and  it  can  be  seen  from  fig¬ 
ures  3  and  4  that  the  expression  for  Inv[G(t))  diverges  as  the  upper  and  lower  integral 
limits  ±a  are  increased  in  wave-number  space.  Nevertheless,  engineering  approxima¬ 
tions  of  /nffG^f)]  are  frequently  used  by  cutoff  of  the  limits  of  the  Fourier  integral 
(also  called  finite-bandwidth  approximations,  since  the  support  of  the  gaussian  filter 
is  reduced  from  ]  —  oo,  +oo[  to  [-a,+a]).  In  other  words,  information  on  the  impulse 
response  and  the  output  for  large  wave-numbers  k  has  to  be  discarded.  However,  the 
approximations  to  the  original  input  thus  obtained  are  highly  sensitive  to  the  cutoff 
value  a ,  and  if  there  is  no  a  priori  knowledge  of  the  Fourier  spectrum  width  of  the 
input,  the  complete  recovery  of  the  original  input  through  Fourier  transforms  is  an 
ill-posed  problem.  Nevertheless,  Gaussians  are  impulse  responses  to  a  wide  variety 
of  physical  systems1,  generally  called  blurring  systems.  Systems  with  pure  Gaussian 
impulse  response  are  responsible  for  gaussian  blur,  while  systems  without  a  Gaus¬ 
sian  impulse  response  are  responsible  for  non-gaussian  blur.  Telescopes  and  a  vast 
majority  of  optical  systems  are  blurring  systems  in  two  spatial  dimensions.  Commu¬ 
nication  satellites  are  blurring  systems  in  one  dimension  (waveforms  are  transmitted 
sequentially).  If  the  convolution  inverse  of  a  Gaussian  (of  infinite  support)  cannot  be 
defined  as  a  smooth  point-function,  how  can  it  be  defined  and  subsequently  utilized 
to  solve  the  inverse  problem? 

We  propose  here  a  novel  approach  to  Gaussian  deconvolution  by  constructing  a  rig¬ 
orous  relation  between  the  deconvolution  of  a  waveform  that  results  from  a  Gaussian 
filter,  and  the  integration  of  the  diffusion  equation  with  negative  diffusi vity  (antidif¬ 
fusion).  In  constructing  this  relationship,  a  new  class  of  highly  singular  functions  that 

'or  at  least  good  approximations  thereof 
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are  called  hyperdistributions  are  introduced,  as  well  as  their  algebraic  properties  and 
a  sequence  of  analyzing  smooth  point-functions  (Hermite-Rodriguez  wavelets)  that 
approximate  them. 


2  Singular  Functions:  Distributions  and  Hyper¬ 
distributions 

2.1  Taylor  and  Moment  Expansions 

The  theory  of  hyperdistributions  is  built  on  ideas  related  to  Dirac’s  delta  function, 
which  is  technically  a  dis/n‘6u<jon[4]-[l  1).  The  Dirac  delta  function,  S(x),  has  the 
following  properties: 

/+oo 

S(x)dx  =  1  (21) 

•oo 

and  Dirac’s  identity: 

f+o o 

/(*)  =  /  S(x-  x')f(x')dx'  (22) 

for  any  suitably  smooth  function  f(x). 

Distributions,  like  the  Dirac  delta  function,  were  first  introduced  by  L.  Schwartz[4], 
and  widely  publicized  by  G.  Temple[6]  and  M.  Lighthill[10].  In  fact,  Lighthill  dedi¬ 
cated  his  book  to  Dirac,  Schwartz,  and  Temple,  in  the  following  way: 

To 

PAUL  DIRAC 

who  saw  that  it  must  be  true, 

LAURENT  SCHWARTZ 
who  proved  it, 

AND 

GEORGE  TEMPLE 

who  showed  how  simple  it  could  be  made 

Both  Dirac  and  Temple  are  Englishmen,  while  Schwartz  is  a  Frenchman.  The  English 
do  have  a  knack  for  putting  down  the  French,  do  they  not? 

Distributions  allow  for  considerable  simplification  and  increased  mathematical  ele¬ 
gance  in  the  handling  of  integral  and  differential  equations.  For  example,  point-forces, 
or  short-time  impulses,  are  frequently  mathematically  modelled  as  Dirac  delta  func¬ 
tions.  All  expressions  involving  distributions,  are  assumed  to  hold  under  integration 
with  any  test  /unctjon(6](10).  More  specifically,  distributions  are  defined  by  a  se¬ 
quence  of  functions,  and  the  property  of  “weak  convergence”.  Good  functions  (total 


functions  that  are  differentiable  to  all  orders,  and  taper  at  ±oo  faster  than  any  power) 
play  the  role  of  “testing”  a  sequence  of  functions  for  weak  convergence.  That  is,  a 
sequence  of  good  functions  {/„(x)}„  is  a  distribution  if 

/+oo 

<f>(x)fn(x)  dx  <  oo  for  all  good  functions  <f>  (23) 

-oo 

For  a  more  detailed  introduction  of  distributions,  refer  to  section  2.4:  ’  Proof  of 
Consistency:  Taylor  Expansion  of  the  Delta  !  unction. 

We  now  expand  S(x'  —  x)  in  a  Taylor  expansion2,  in  terms  of  its  derivatives  about 
x': 

6(x'  -x)  =  S(xf)  -  xS\x')  +  ~S’’(x')  +  ...  (24, 

where  6"(x)  =  ^ Six).  A  “local”  approximation  to  fix)  can  be  derived  from  Eq.(24) 
in  the  following  way.  By  substituting  Eq.(24)  into  Dirac  ’s  identity,  we  recover  the 
usual  Taylor  expansion  of  /(x)  about  x  =  0: 

fix)  =  m  +  xf'i  0)  +  |^/"(0)  +  ...  +  Rn{x)  (25) 

T  his  approximation  is  local  in  the  sense  that  it  requires  derivatives  of  fix)  at  a  single 
point,  x  =  0,  and  in  general  has  a  limited  radius  of  convergence.  On  the  other  hand, 
by  expanding  6(x  —  x')  in  a  Taylor  expansion  about  x: 

<5(x  -  x')  =  6(x)  -  x'S'ix)  +  —  <5”(x)  +  ...  (26) 

When  this  series  is  substituted  into  Dirac ’s  identity,  We  obtain 

A/2 

fix)  =  M°S(x)  -  A/V(x)  +  -fi-S  '(*)  +  T  Rnix)  (27) 

The  coefficients  A/n,  defined  by 


/+oo 

x”/(x)dx  (28) 

•OO 

are  the  centered  moments  of  the  function  fix).  Therefore,  Eq.(27)  is  an  approxima¬ 
tion  of  fix)  involving  global  information  about  fix). 

The  global  expansion  of  /(x)  in  terms  of  its  centered  moments  and  derivatives  of 
the  Dirac  delta  function  is  the  motivation  for  the  introduction  of  hyperdistributions. 

2 1 li is  is  accomplished  in  a  formal  way;  the  justification  and  proof-of-consistencv  will  he  given  in 
section  2.4:  *  Proof  of  Consistency:  Taylor  Expansion  of  the  Delta  Function 
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Definition  (formal):  Any  “function”  which  may  be  formally  written  a s, 

/(*)  =  f>„«"(x)  (29) 

n=0 

with  finite  and  real  coefficients  an  given  by 

an  =  (— l)n— r-  with  Mn=  [+°°  xnf(x)dx  (30) 

Tt\  J — oo 

defines  a  hyperdistribution. 

To  carry  on  the  analogy  with  natural  numbers  introduced  previously,  hyperdis¬ 
tributions  can  be  identified  with  the  set  of  negative  integers.  Much  in  the  same  way 
hyperdistributions  are  not  smooth  point-functions,  negative  integers  do  not  corre¬ 
spond  to  any  physical  reality.  Indeed,  we  can  physically  represent  the  integer  2  by 
picturing,  say,  two  apples  (or  physical  entities  of  any  other  sort)  in  our  mind.  But 
how  do  we  represent  the  negative  integer  -2?  Is  an  “absence”  a  constructive  defini¬ 
tion?  In  spite  of  this  conceptual  difficulty,  negative  integers  are  required  in  simple 
arithmetic,  since  they  are  inverses  to  the  set  of  natural  numbers  for  the  operation 
of  addition.  That  is  the  reason  why  the  space  of  natural  numbers  is  extended  to 
the  space  of  positive  and  negative  integers.  In  a  similar  fashion,  hyperdistributions 
do  not  correspond  to  any  physical  reality  (e.g.  hyperdistributions  are  not  smooth 
-or  even  piecewise  smooth-  point-functions),  yet  it  will  be  shown  that  hyperdistri¬ 
butions  are  convolution  inverses  to  smooth  point-functions,  and  are  thus  required  in 
the  “convolution  calculus”  of  signals  and  systems.  One  is  thus  naturally  led  to  ex¬ 
tend  the  space  of  smooth  point-functions  to  include  highly  singular  functions  such  as 
hyperdistributions. 


2.2  The  Convolution  Group 

Given  any  two  hyperdistributions  f\  and  /2,  their  linear  combination  A/i  +  /t /2  is  also 
a  hyperdistribution  (where  A  and  p  are  real  numbers): 


/i(x)  = 
/*(*)  = 
A/i(x)  +  m/2(x)  = 

The  pth  derivative  dpf(x)/dxp  =  Vp/(x) 
distribution: 

OO  OO 

V7(x)  =  Vp  £  an6n(x)  =  Y,  M”(X 


OO 


Y  a^n(x) 

(31) 

n= 0 

oo 

£  W(x) 

(32) 

n=0 

oo 

Y(^an  +  pbn)Sn{x) 

(33) 

n=0 

f  a  hyperdistribution  f(x)  is  also  a 

hyper- 

where  6=1°  if  "  <  P 

„  [  a„.p  otherwise 

(34) 
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The  convolution  of  two  hyperdistributions  f\  *  f2  is  a  hyperdistribution: 


/!(*)•/»(*)  =  I>„V"*(l)  .  £>„V"«(l) 


(35) 


«=o 


n=0 


-  r (E'hP**W)(Z. >.v;_,.«(x-x')),iz'  (36) 

■J—GO _ r* 


P=0  9=0 

Noting  that  VJ_X,  =  (  —  1)’VX,,  and  with  9  integrations  by  parts,  we  obtain: 

00  00  /-f  OO 

/iW  *  M*)  =  (-J)1’  £(«»£  5,  /  V?’6(x')6(x  -  x')dx') 

p=0  p=o 

From  Dirac’s  identity, 

/+°°  V^96(x')6{x  -  x')dx'  =  Vp+,<5(*) 
and  thus  Eq.(36)  may  finally  be  reduced  to 


(37) 


(38) 


OO  OO 


«*)  *  A(x)  =  EKL^V’^x)) 

p=0  ?=0 


(39) 


or,  equivalently,  with  r  =  p  -fi  q: 


/.w*/*(*)=e(e^  r^) 

r=0  fp=0  J 


(40) 


Since  {a„}„  and  {6„}n  are  real  and  finite,  {XTp=0  ap^r-p}  are  real  and  finite,  and 
Eq.(40)  defines  a  hyperdistribution.  The  three  properties  in  Eqs.(33),  (34),  and  (40) 
above  may  be  thought  of  as  “closure  properties”  of  hyperdistributions. 

Another  attractive  property  of  hyperdistributions  is  their  behavior  under  the 
Fourier  transform.  Again  taking  f(x)  =  GnVn<5(x),  the  Fourier  transform  can 
be  evaluated  as: 


/+°°  00  r+00 

e-'kxf{x)dx  =  X>"/  e~'kxVn6(x)dx 

■OO  _ n  J  —  OO 


n=0 

00 


=  E  “«(-!)’ /  “(V”c-“*)«(x)dr 

OO 

=  £«.(*’*)" 


(41) 

(42) 

(43) 


n=0 


where  we  have  used  integration  by  parts  n  times,  and  Dirac’s  identity  with  Vnc,kx  |r-0  = 
(ik)n.  The  Fourier  transform  of  a  hyperdistribution  is  thus  a  formal  power  series  in 
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the  wave-number  k.  Moreover,  the  Fourier  transform  of  a  function  is  the  “moment 
generating  function”,  because: 


A/(n)  ,  ,  r°° 

an  =  (-!)" — —  with  M{n)=  xnf(x)dx  (44) 

72:  J — oo 

As  a  result,  the  definition  of  hyperdistributions,  which  requires  that  the  coefficients 
an  be  real  and  finite,  is  equivalent  to  requiring  that  its  Fourier  transform  be  a  real 
analytic  function  of  the  wave-number  k. 

Hyperdistributions  allow  an  effective  computation  of  the  convolution  inverse.  Given 
a  hyperdistribution  /,  the  desired  convolution  inverse  Inv[f]  satisfies 

f*Inv[f]  =  S  (45) 

where  S  represents  Dirac  ’s  delta  function,  which,  by  Dirac’s  identity,  is  the  unit 
element  of  the  convolution  operation.  We  shall  show  by  construction  that  if  /  can  be 
written  as  a  hyperdistribution,  then  Inv[f]  is  also  a  hyperdistribution.  Bv  definition, 

(f  *  Inv[f])(x)  =  f  f(x')Inv[f](x  -  x')dx' =  6(x)  (46) 

It  is  now  necessary  to  compute  the  product  of  the  sums  and  match  coefficients.  Taking 
/(x)  —  5Z„anV"6(x),  and  Inv[f(x)]  =  £n  6„Vn£(x),  it  is  seen  that  the  computation 
of  the  convolution  inverse  is  effectively  the  determination  of  a  collection  of  bn  values, 
given  a  set  of  an  values.  Substituting  into  Eq.(46), 

/+oo  °°  °° 

(£  ^(x'))(£  bgVqS(x  -  x'))  dx'  =  S(x)  (47) 

p= 0  1=0 

Once  again  noting  that  V’_x,  =  (  —  lJ’Vx-  and  that  Vp6(x)  *  Vq6(x)  =  Vp+,<5(x), 
Eq.(47)  may  finally  be  reduced  to: 


E 


r= 0 


Vr<5(x)  =  6(x) 


(48) 


Matching  coefficients  on  the  left-hand  and  right-hand  sides  implies  that  only  the  r  =  0 
term  survives.  The  result  is  a  linear  system  of  equations  for  the  b  values  in  terms  of 
the  a  values.  It  is  easiest  to  see  the  behavior  by  writing  the  first  few  equations  in  this 
linear  system, 


ao&o  =  1 

a0bi  +  Gj  bo  =  0 

a0^2  +  01^1  -f  G260  =  0 

G0&3  +  d\b2  +  02^1  +  03^0  =  0  (49) 
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and  so  forth.  Thus,  it  can  be  seen  that,  bo  =  1/ao,  b i  =  —a^/a^  b2  =  cl\/(1q  —  a2/al, 
etc..  The  linear  system  in  Eq.(49)  is  called  the  Bochner  algebra  for  hyperdistributions, 
after  a  very  similar  algebra,  developed  by  S.  Bochner[12]  for  power  series. 


In  matrix  notation,  and  if  a0  =  1,  the  coefficients  bn  are  given  by  the  following 
Toeplitz[13]  determinant: 


bn  =  (-I)” 


a  i 
a2 

“3 


1 

ai 

a2 


0 

1 

ai 


0 

0 

1 


Q-n  &n  —  1  ^n-2  3 


0 

0 

0 


1 

«1 


(50) 


We  note  that  Inv[f  ]  is  indeed  a  hyperdistribution,  and  is  determined  by  this  algorithm 
up  to  a  function  with  vanishing  moments. 

Some  convolution  inverses  appear  in  the  familiar  context  of  potential  theory.  For 
example,  the  "isotropic  quadrupole”,  which  is  a  distribution,  is  the  convolution  inverse 
of  the  Coulomb  potential  in  3  dimensions: 

V2<S(r  )  =  Inv[—  ] 

4xr 

Other  simple  examples  of  convolution  inverse  pairs  in  one  dimensional  potential  theory 
are  (|  x  |,<5"(i)),  (sgn(x),6'(x)),  and  (6(i),6(x)).  More  examples  of  convolution 
inverses  can  be  obtained  from  standard  Green’s  functions.  For  example,  the  one- 
dimensional  Helmholtz  equation 

<P 

u  +  p2u  =  0  (51) 

dxz 

can  be  analyzed  with  hyperdistributions,  as  well  as  with  Fourier  transforms,  to  obtain 
both  its  Green’s  function  and  the  convolution  inverse  of  its  Green’s  function.  Thus, 

(i?+i)e'w  =  ~4,l)  <52) 

The  Bochner  algebra  for  hyperdistributions  shows  that  only  the  terms  b0  and  b2 
contribute  to  the  convolution  inverse  of  the  Green’s  function,  thus  yielding: 

Inv[e-'*'/2}  =  ^(z)  -  6"(x)  (53) 

The  result  above  is  easily  reprodu^d  by  utilizing  Fourier  transforms.  However,  it  was 
shown  in  Eq.(19)  that  the  convolution  inverse  of  a  Gaussian  cannot  be  obtained  with 
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Fourier  transforms.  In  contrast,  that  same  convolution  inverse  can  easily  be  obtained 
through  the  Bochner  algebra  for  hyperdistributions:  calculation  of  the  power  moments 
of  a  Gaussian  £a(z)  of  width  A 

<a(i)  =  (54) 

yields: 

/+OO 

x2ne  'Z/yfnXdx 

■OO 

_  A2,  (2n-l)V. 

2” 

a2n+l  =  0 

where  the  double  factorial  involves  exclusively  the  product  of  odd  numbers.  This 
leads  to  the  hyperdistribution  expansion  of  6\(x): 

00  /  A  \2n  1 

-  s  y  <57) 


The  Bochner  algebra  for  hyperdistributions  yields  the  convolution  inverse  of  the  Gaus¬ 
sian  6\{x): 


OO 

Inv[Sx{  x)]  = 

71=0 

(y 

(58) 

In  a  more  compact  operator  notation: 

S\[x) 

(591 

Inv[6\(x)} 

=  S(t) 

(60) 

In  this  section,  it  was  shown  that  the  convolution  of  any  two  hyperdistributions 
yields  another  hyperdistribution.  Furthermore,  it  was  proven  that  any  hyperdistri¬ 
bution  has  another  hyperdistribution  as  its  convolution  inverse,  provided  the  impulse 
response  can  also  be  written  as  a  hyperdistribution.  A  general  algorithm  for  decon¬ 
volving  any  hyperdistribution  -the  Bochner  algebra-  was  given.  As  a  result,  hyper- 
distributions  provide  a  closure  for  the  classical  semi-group  of  smooth  point-functions 
with  the  operation  of  convolution[14|.  If  7i  denotes  the  space  of  Hyperdistributions, 
with  the  space  of  smooth  point  functions  embedded  in  it,  (7-f,*)  is  an  abelian  group, 
with  the  Dirac  delta  as  the  unit  element.  Hyperdistributions  are  thus  an  extremely 
useful  tool  for  solving  integral  equations  of  convolution-type[15|. 


(55) 

(56) 
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2.3  Function-Theoretic  Definition  of  Hyperdistributions 


The  process  for  defining  hyperdistributions  parallels  the  Temple  definition3  of  a  dis¬ 
tribution  (generalized  function)  as  a  good  sequence  of  good  functions[4]-[llj.  Good 
functions  are  smooth  and  tapered.  More  precisely,  they  are  total  point  functions  that 
are  differentiable  to  all  orders  (C°°),  and  decay  at  ±oo  faster  than  any  power.  Good 
functions  play  the  role  of  “testing”  a  sequence  of  good  functions  for  weak  convergence. 
In  fact,  a  sequence  of  good  functions  is  a  distribution  if 

/+oo 

4>(x)fn(x)dx  <  oo  for  all  good  functions  <f>  (61) 

•oo 

Since  hyperdistributions  are  conceived  as  “generalized”  distributions,  a  second  order 
generalization  of  functions  is  in  fact  implemented.  Consequently,  a  double  test  is 
needed  as  a  convergence  criterion.  This  criterion  is  implemented  by  introducing  very 
good  functions  Ga(^)  with  the  following  properties: 

1.  G\(x)  is  smooth,  that  is,  differentiable  to  all  orders,  i.e.  C°° . 

2.  G\(x)  is  essentially  compact,  i.e.  it  has  a  Gaussian  decay  at  ±oo: 

Ga(z)  ~  limx_00  Ne~x*ls2 . 


12 


The  sequence  {'W^A.n,  where  A  is  a  nonnegative  real  and  n  is  a  natural  number,  is  a 
good  sequence  if,  for  all  good  functions  <f>  and  for  all  very  good  functions  Ga,  there 
exists  a  A0  such  that,  for  all  A  >  Ao, 


/+°° 

^(x)(^n  *  G\){x)dx  <  oo  (66) 

A^U  00 

We  note  that  e±tVl S(x  )  are  hyperdistributions.  The  sum  (hyperdistribution) 

(67) 

k= 0 

can  thus  be  viewed  either  as  a  sequence  of  “good”  distributions  as  n  — >  00 

t.akVkS(x)  (68) 

fc= 0 

or,  as  A  — ►  0,  as  a  sequence  of  good  functions: 

OO 

£a*V%(x)  (69) 

Ar=0 


The  latter  representation  is  a  Rodrigue:  expansion.  The  Rodriguez  formula  for  Iler- 
mite  polynomials[2]  can  be  used  to  show  that  the  derivatives  of  a  Gaussian  form 
a  complete  set  of  orthogonal  polynomials  in  an  L 2  space.  And  thus  the  Rodriguez 
expansion  yields  a  very  useful  point  function  approximation  to  any  hyperdistribution: 

00  r  j  /  x  \ 

<70> 

k=o  * 

where  //„(x)  denotes  the  Hermite  Polynomial  in  x  of  order  n. 

The  problem  of  “approximating”  hyperdistributions  with  smooth  point-functions 
will  be  studied  shortly. 

2.4  *  Proof  of  Consistency:  Taylor  Expansion  of  the  Delta 

Function 

We  start  by  quoting  three  useful  definitions  utilized  by  LighthillflO]  in  his  definition 
of  distributions4. 

4Lighthi!i  actually  uses  the  terminology  “generalised  function",  instead  of  “distributions",  which 
was  first  introduced  by  Schwartz[4] 


13 


Definition  1  (p.  15).  A  good  function  is  one  which  is  everywhere  differ¬ 
entiable  any  number  of  times  and  such  that  it  and  all  its  derivatives  are 
0(1  x  |-JV)  as  |  x  |— f  oo  for  all  N. 

Example  1.  e~x  is  a  good  function. 

Definition  2.  A  fairly  good  function  is  one  which  is  everywhere  differen¬ 
tiable  any  number  of  times  and  such  that  it  and  all  its  derivatives  are 
0(|  x  |-iV')  as  |  x  |— *  oo  for  some  N. 

Example  2.  Any  polynomial  is  a  fairly  good  function. 

Definition  3.  A  sequence  fn(x)  of  good  functions  is  called  regular  if,  for 
any  good  function  F(x)  whatever,  the  limit 

/+oo 

fn(x)F(x)dx  (71) 

•oo 


exists. 

(...)  Definition  5.  A  generalised  function  f(x)  is  defined  as  a  regular 
sequence  f„(x)  of  good  functions  (...).  The  integral 


f(x)F(x)dx 


72) 


of  the  product  of  a  generalised  function  f(x)  and  a  good  function  F{ x)  is 
defined  as 


fn(x)F(x)dx 


(73) 


In  other  words,  Lighthill  defines  distributions  by  a  sequence  of  good  functions 
(total  functions,  differentiable  to  all  orders,  that  taper  at  ±oc  faster  than  any  power), 
and  by  the  property  of  weak  convergence,  which  requires  that  the  limit  in  Eq.(71) 
exists,  for  any  choice  of  good  (testing)  function  F(x). 

We  now  justify  the  expansion  in  Eq.(  24)  of  S(x'  —  x)  in  a  Taylor  series,  in  terms 
of  its  derivatives  about  x'  (since  6(x'  —  x)  is  technically  a  distribution,  and  not  a 
point-function,  such  a  justification  and  proof-of-consistency  is  required).  The  Taylor 
expansion- with-remainder  of  S( x'  —  x)  about  x': 


,,v-i 


Hr'  -r)  =  Hr')  +  xi'(r')  4  £<V)  4  ...  +  +  RS  (T1 


7’o  be  consistent  with  Lighthill’s  definition  of  distributions,  We  are  required  to  test 
Kq.(  71)  for  aj]  available  good  functions  d>(x): 


/  +  <■> o  /■+o O  f  t  J-n  ,  „) 

6{x'  -  x)d(x')dx'  —  J  |  y:  —  S(n,(x')  +  t{s  |  &(x')d.r' 


(75) 
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By  Dirac’s  identity  (Eq.(  22)), 


1  r  -f-OO 

<f>(x)  =  5Z  — r<^(n)(0)  4-  nN  where  KN  =  /  RN <f>{x')dx'  (76) 

n=0  n!  -7-00 

and  we  recognize  Taylor’s  expansion  for  the  test  function  <£(x).  If  such  an  expansion  is 
indeed  legitimate  for  all  test  functions  <t>(x),  then  Eq.(  74)  is  justified  and  consistent. 
However,  Lighthill  does  not  require  of  test  functions  that  their  Taylor  expansions 
converge,  or,  in  other  words,  that  all  test  functions  be  real-analytic. 

Theorem  of  Consistency  for  Hyperdistributions :  In  restricting  the  space  of 
test  functions  to  good  functions  that  are  real-analytic  (i.e.  C"),  Eq.(  74) 
is  justified,  and  the  introduction  of  the  hyperdistribution  sum  in  Eq.(  29) 
is  consistent. 


It  might  seem  paradoxical  that  by  restricting  the  space  of  Lighthill  test  function  (from 
C°°  to  Cw),  we  expand  Lighthill ’s  space  of  singular  functions  to  include  hyperdistribu¬ 
tions,  which  are  even  more  singular  than  distributions.  In  fact,  this  is  a  consequence 
of  the  “dual”  behavior  of  Lighthill’s  regular  sequences  and  Lighthill's  test  functions 
under  the  weak  convergence  property.  Since  the  limit 


fn(x)F{x)dx 


(77; 


must  exist,  extending  the  space  of  regular  sequences  fn{x )  necessarily  restricts  the 
space  of  available  test  functions  F(x). 


2.5  Approximating  Hyperdistributions  with  Smooth  Point- 
Functions 

In  this  section,  we  will  parallel  Lighthill’s  approach  in  defining  the  delta  function  by 
a  sequence  of  narrowing  gaussians.  We  will  introduce  a  sequence  of  smooth  point- 
functions,  whose  limit  is  in  fact  a  hyperdistribution.  These  smooth  point-functions 
will  then  be  the  focus  of  section  3:  Hermite-Rodrigucz  Wavelet  Analysis. 


2.5.1  “Diffusing”  Hyperdistributions 


It  is  observed  that  by  formally  convolving  a  hyperdistribution  flnWA(r)  with  a 
Gaussian  of  width  A: 


<M*) 


(78) 
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one  obtains  an  infinite  series  of  derivatives  of  Gaussians: 


oo  oo  oo 

E  a„V"*(x)  *  6x(z r)  =  £  G«{V"<5(*)  *  *a(*)}  =  E  ^  V^A(x)  (79) 

n— 0  n=0  n=0 

The  first  equality  is  a  formal  one.  The  proof  of  the  second  equality  follows  from  the 
application  of  the  chain  rule  to  the  operator  n  integration  by  parts,  and  Dirac’s 

identity: 


V*6(x)  *  6x(x) 


/+oo 

Vnx_x.6(x-x')6x(x')dx' 

“OO 

(_l)n  f+0OVnxlS(x-x')6x(xf)dx' 

J  —  OO 
_  r+oo 

(-l)2n  /  6(x-x')V$6x(x’)dx' 

J — OO 

Vn6x(x) 


(80) 


In  other  words,  diffusion  “maps”  hyperdistributions,  which  are  singular  functions, 
into  series  of  smooth  point-functions.  This  result  is  not  surprising.  Diffusion  is  a 
smoothing  operation,  and  is  expected  to  map  singular  functions  into  smooth  functions. 
This  is  indeed  the  case  for  Schwartz  distributions:  Picture  a  single  Dirac  delta  as  the 
initial  condition  (at  t  =  0)  of  a  diffusion  process  in  time.  The  1  +  1  homogeneous 
diffusion  equation  can  be  written  as 

|-V’u  =  0  (81) 

The  initial-condition  Green’s  function  G*(x)  for  such  a  diffusion  equation  is  the 
Gaussian[22]  G((x)  =  S2yj(x).  In  other  words,  if  f(x)  is  the  initial  condition  to 
the  diffusion  process,  the  solution  of  the  diffusion  equation  at  time  t  can  be  written 

as: 

u{x,t)  =  Gt{x)  *  J{x)  =  82ft  *  /(x)  (82) 

or,  in  operator  form [23], 

u(x,!)  =  e‘v  ’/(x)  (83) 

Since  in  this  case  f{x)  =  6(x),  we  have 

u(x,<)  =  Gt(x)  *  6(x)  =  S2yi(x)  *  6(x)  =  S2y-e(x)  (84) 

The  last  equality  follows  once  again  from  Dirac’s  identity.  In  other  words,  at  i  =  c, 
for  any  <  >  0  however  small,  a  Dirac  delta  is  “mapped”  into  a  Gaussian  of  width 
VTt  by  the  diffusion  equation.  Similarly,  it  can  be  proven  using  Eq.(80)  that  a  finite 
derivative  of  a  Dirac  delta,  gets  mapped  to  the  derivative  of  same  order  of  a  Gaussian 
of  width  \fTl,  by  the  diffusion  equation: 

u(x,t)  =  Gt(x)  *  VpS(x)  -  S2yi(x)  *  V”S(x)  =  Vps2yt(x)  (85) 
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We  now  return  to  hyperdistributions,  by  picturing  an  infinite  number  of  Dirac  delta 
functions,  all  centered  at  the  origin,  and  with  different  amplitudes  a„,  that  is:  an^7n<5(x), 

as  an  initial  condition  to  a  diffusion  process.  That  hyperdistribution  will  formally  dif¬ 
fuse  at  time  t  >  0,  however  small  t  may  be,  into  the  following  functional  series: 

f (86) 

n=0 

In  other  words,  diffusion  maps  hyf  erdistributions  into  series  of  smooth  point  functions. 

If  the  series  converge,  we  can  say  that  diffusion  maps  hyperdistributions  onto  the  space 
of  smooth  point-functions.  The  sufficient  conditions  f<~"  convergence  of  the  series  in 
Eq.(86)  remains  however  to  be  investigated. 

To  this  end,  the  Rodriguez  formula  for  Hermite  Polynomials[2]: 

(-ir//„(z)e-"2  =  Vne-r2  (87) 

can  be  rescaled  as  follows: 

=  V"6,(i)  (88) 

where 

«.*(*)  =  ^  (89) 

And  thus  Eq.(86)  becomes  a  series  involving  Hermite  Polynomials: 

OO 

£>„(— i  rfrtiiy.)  (%) 

n=0 

Unfortunately,  the  series  above  is  not  a  classical  Hermite  series,  since  the  coefficients 
a„  are  centered  power-moments  given  by  Eq.(  30)  (weighted  by  monomials  in  x, 
and  not  by  Hermite  polynomial)  As  a  result,  one  cannot  use  the  standard  tools 
of  Christoffel-Darboux  theory[2],  which  yield  the  sufficient  conditions  for  the  conver¬ 
gence  of  Hermite  series[2].  The  task  that  lies  ahead  is  to  transform  t>  e  series  in  Eq.(90) 
into  standard  Hermite  series  by  relating  the  centered  power- moments  ffff  xnf(x)dx 
to  the  centered  Hermite  moments  ffff  Hn(x)f(x)dx,  where  Hn{x)  denotes  the  Her¬ 
mite  polynomial  of  order  n.  That  is  the  purpose  of  the  C-matrix  transform,  which 
will  be  studied  shortly. 

2.5.2  Hermite-Rodriguez  Expansions 

We  now  expand  the  function  f(x)  directly  in  terms  of  Hermite  polynomials.  Hermite 
polynomials  constitute  an  orthogonal  set  of  basis  functions  for  the  appropriate  L2 
space,  and  such  an  expansion  is  thus  legitimate.  We  start  with  the  expansion  of  f(x) 
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in  terms  of  derivatives  of  the  Gaussian  of  width  A,  with  yet  unknown  coefficients  (or 
moments)  bn: 

/(*)  =  £  KV"h(x)  (91) 

n=0 

where 


6x(x)  = 


v/tFA 


(92) 


With  the  use  of  the  Rodriguez  formula  for  scaled  Hermite  polynomials  derived  in 
Eq.(8S): 

(-l)”//„A(i)«x(x)  =  V"^(i)  (93) 

with 

//.*(*)  =  (94) 

We  substituteVn<5A(2:)  in  Eq.(91)  with  the  expression  in  Eq.(93).  Eq.(91)  becomes: 


+oo 


/<x)  =  x:M-ir//n1(x)Wx) 


n= 0 


We  now  utilize  the  orthogonality  of  Hermite  polynomials: 

[+°°  Hn(x)Hm(x)e-l7dx  =  yft2nn\8n 

J  —  oo 


(95) 


(96) 


where  <5nm  represents  the  Kronecker  delta,  which  is  0  if  n  ^  m  and  1  if  n  =  m. 
We  multiply  Eq.(95)  with  H^(x)  and  integrate  from  -oo  to  +oo.  Assuming  the  sum 
in  Eq.(95)  converges,  We  can  interchange  the  integral  and  sum  symbols,  and  utilize 
Eq.(96)  to  obtain: 

(-l)X  =  /+“ /(x)//n*(x)dx  (97) 

Z  Tl.  J — oo 

The  coefficients  bn  above  are  called  the  Hermite- Rodriguez  moments  of  the  Hermite- 
Rodriguez  expansion  of  f{x): 


/(*)  =  D-iri,«,,(I)(*) 

n~Q 


(98) 


hereafter  denoted  as  the  HR  expansion  of  j(x)  with  weight  A.  The  weight  parameter 
is  a  novel  feature  of  HR  expansions,  when  compared  to  standard  expansions  in  terms 
of  complete  sets  of  basis  functions.  The  standard  expansions  do  not  contain  a  free 
parameter.  This  parameter  may  be  employed  to  speed-up  the  convergence  of  HR 
expansions.  In  other  words,  for  every  function  f(x)  to  be  approximated  by  an  HR 
expansion,  there  exists  one  (or  possibly  more)  choice(s)  of  the  weight  A  for  which  the 
HR  expansion  converges  to  a  good  approximation  of  f(x)  with  a  minimum  number 
of  partial  sums. 


2.5.3  Power  Moments,  Hermite-Rodriguez  Moments,  and  the  C-Matrix 
Transform 


We  now  relate  the  centered  power  moments  of  a  function  f(x): 

(-l)n  r+oo 
_  _  V  >  I  _n  rt 


—  1  r+oo 
fi !  J — oo 


with  the  Hermite-Rodriguez  moments  of  the  same  function: 


.»  (-l)nA3n  /•+» 


‘A  r+°°  > 

-j—  f  (x)f(x)dx 

Til  J  —  oo 


(100) 


By  utilizing  the  classical  power-series  expansion  of  Hermite  polynomials[2]: 

ln/2]  nn—2p 

^,=n!5(-1)p^)T^ 

where  [n/2]  denotes  the  natural  number  inferior  or  equal  to  the  rational  n/2, 
substituting  Eq.(lOl)  in  Eq.(lOO),  we  obtain: 


/  I  \n  \n  [«/2] 


2"— 2p  r+oo 


p}(n  —  2p)!  J- oo 


/+oo 

(z/A)"~2p/(z)  dx 

-oo 


and  thus, 


[n/2]  \2p 


(101) 


and  by 


(102) 


(103) 


In  matrix  notation: 


*f  _ 

k ^  —  C  . 


where  the  matrix  Cx  is  lower-triangular,  and  is  given  by: 

__  /  (-1)^  if  (*  -  j)  is  even  and  positive 

3  1  0  otherwise 


(104) 


(105) 


The  determinant  of  Cx  is  always  equal  to  1.  As  a  result,  one  can  always  transform  an 
HR  expansion  of  /(z)  into  a  power-moment  series  expansion  of  the  same  function,  and 
vice-versa.  The  problem  of  convergence  of  the  series  in  Eq.(90)  has  effectively  been 
reduced  to  the  problem  of  convergence  of  Hermite  expansions.  ChristofTel-Darboux 
thcory[2]  then  yields  sufficient  conditions  for  convergence  of  Hermite-Rodriguez  ex¬ 
pansions. 
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2.5.4  Approximating  Hyperdistributions  with  Hermite-Rodriguez  Expan¬ 


sions 


/(x)  is  expanded  as  a  hyperdistribution: 


/(*)  =  i) 


(106) 


where 


(  — l)n  r+°° 

an  =  ~~  *V(*) 

n !  J —oo 


(107) 


That  hyperdistribution  is  then  approximated  with  a  Hermite-Rodriguez  series  of 
weight  A: 

/M  =  i 'XV'M*)  (108) 

n=0 

where  the  Hermite-Rodriguez  moments  bn  are  given  by  the  two  equivalent  equations 
below: 


bt  =  5-^^/+”/(x)Wn'(x)rfx 
f-l)n  /’+°° 

«i  =  /  **/(*) * 

n!  j — oo 

with  6*  =  (C*m)  am 


K  = 


(109) 


(110) 


The  sufficient  condition  for  convergence  of  Hermite-Rodriguez  expansions  is  given 
in  the  following  section. 


2.5.5  Convergence  of  Hermite-Rodriguez  Expansions 

ChristofTel-Darboux  theory(2]  yields  the  sufficient  condition  for  convergence  of  Her- 
mite  series: 

°°  1  f+°°  7 

£**.(*)  with  Cn  =  -~—=  e~x  f(x)IIn(x)dx  (111) 

n=o  ^  nVrr  J- oo 


converges  to  /(x)  if  the  following  condition  is  met[2): 


/  +  °°  7  „ 

e~x  j2(x)dx  <  oo 

-OO 


(112) 


Hermite-Rodriguez  expansions  arc  of  the  form 


I'n(x/\) 


\2n  r+oo  lL(x/\ ) 

6x(x)  with  bn  =  ±-f  f{x)L^ldx  (113) 

/  n!  J  —  oo  A 
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By  multiplying  Eq.(lll)  by  and  by  rescaling  from  x  to  x/A,  the  convergence 
criterion  for  Hermite  Rodriguez  expansions  is  obtained.  It  is  the  following  sufficient 
condition.  The  series  in  Eq.(113)  converges  if  the  following  condition  is  met: 

f  erVA2y2(x)  dx  <  oo  (114) 

J  —  OO 

2.6  Hyperdistributions  in  2  Dimensions 

The  generalization  of  hyperdistributions  to  two  spatial  dimensions  is  now  introduced. 
Formal  expressions  of  the  type 

OO  OO 

/(*.»)  =  El ;  a™  v;v™  «(*,„)  (115) 

Ti— 0  m—Q 

which  are  called  nyperdistributions,  form  an  algebraic  field.  V”  and  V”  respectively 
denote  the  nth  and  mtK  partial  derivative  operators  in  the  x  and  y  directions,  and 
6(x,y)  denotes  the  2D  Dirac  delta  function: 

6(x,y)  =  S(x)6(y)  (116) 

The  coefficients  anm  are  given  by 

(-  rr+oo 

dnm  =  - j — f—  //  xnymf(x,y)dxdy  (117) 

72.771.  JJ — oo 

If  /(x,y)  is  convolved  with 

OO  OO 

^y)  =  EE6-V"V^«5(x,y)  (118) 

n=0  m=0 

through  the  usual  convolution  product  *,  the  result  is  once  again  a  hyperdistribution, 
with  coefficients  given  by 

c  =  a*  b  (119) 

This  is  the  discrete  convolution  product,  which  can  be  written  explicitly  as 

n  m 

^nm  —  n— p ,  m—qbpq  (120) 

p=0  q= 0 

To  prove  the  above  formula,  we  observe  that 

(V"V^(x,y))  *  (V£V’*(x,y))  =  V"+p  V-+’<5(x,y)  (121) 

which  can  be  verified  by  repeated  integration  by  parts.  The  algebra  represented  by 
Equation  (120)  is  the  2D  Bochncr  algebra\  12].  This  constitutes  an  efficient  means  to 
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deconvolve  2D  images.  In  particular,  the  coefficients  of  the  convolution  product  are 
given  by 


0OO&OO  = 

Coo 

aioboo  4-  floo^io  = 

C\o 

fl0lf*00  +  aoof>oi  = 

C01 

floo^n  +  aoi^io  +  aiof>oi  +  an^oo  — 

Cll 

a20^00  4"  GlO^lO  4  U 00^20  = 

C20 

a02^oo  +  Ooiboi  +  aoobo2  = 

C02 

... 

(122) 

The  algebra  is  straightforward,  as  one  notices  that  each  ctJ 

is  the  sum  of  all 

configu- 

rations  in  ( p ,  q,r,s)  of  the  products  apqbTa,  where  p  -+■  r  =  i  and  q  +  s  =  j. 

If  a  deconvolution  is  to  be  carried  out  (i.e.  to  determine  the  original  image  after 
it  was  convolved  with  a  filter  function),  the  coefficients  bXJ  and  c,j  are  known.  Indeed, 
the  coefficients  are  respectively  the  2D  moments  of  the  filter  function  and  those  of  the 
degraded  image.  This  is  a  straightforward  generalization  to  two  independent  variables 
of  similar  results  in  one  dimension.  The  unknown  coefficients  atJ  are  determined  by 
the  linear  system  in  Eq.(122).  Their  solution  can  be  obtained  rapidly  with  Gaussian 
elimination.  Once  the  moments  a,j  of  the  original  image  are  determined,  the  image 
can  be  reconstructed  by  approximating  the  resulting  hyperdistribution  with  an  IIR 
expansion,  much  in  the  same  way  it  is  done  in  one  dimension.  If,  on  the  other  hand, 
one  is  simply  looking  for  the  convolution  inverse  to  /,  that  is  to  determine  <7  such  that 
f  *  g  =  6(*,y),  then  coo  =  I  and  Cnm  =  0  for  all  (n,m)  ^  (0,0);  and  the  coefficients 
btJ  are  once  again  easily  determined. 

3  Hermite-Rodriguez  Wavelet  Analysis 

3.1  Introduction 

Wavelets  are  an  exciting  new  technique  for  solving  difficult  problems  in  mathemat¬ 
ics,  physics,  and  engineering.  Applications  are  as  diverse  as  seismic  exploration, 
data  compression  of  digitized  signals  and  images,  the  detection  of  submarines,  and 
improvements  in  CAT  scans  and  other  medical  imaging  technologies.  Wavelets  al¬ 
low  complex  information  such  as  speech,  music,  photographs,  or  video  images,  to 
be  broken  down  into  fundamental  building  blocks  -the  wavelets-,  and  subsequently 
reconstructed  with  high  precision. 

A  recent  article  in  Business  Week  magazine  (Feb.  3  1992),  entitled  ‘Wavelets  arc 
causing  ripples  everywhere  ’,  emphasizes: 
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Catching  the  Wavelet-  A  dazzling  theory  of  mathematics  is  leading  to  ad¬ 
vances  in:  data  compression  (Wavelets  identify  key  features  of  an  image, 
allowing  engineers  to  reconstruct  a  photo  with  only  a  tiny  fraction  of  the 
information  in  the  original),  scientific  calculation  (Wavelets  may  help  de¬ 
cipher  phenomena  such  as  turbulence-  leading,  for  example,  to  improved 
aircraft  designs),  signal  analysis  (the  new  method  makes  it  easier  for  mili¬ 
tary  radar  to  spot  hidden  tanks),  and  medical  imaging  (scientists  are  using 
Wavelets  to  produce  magnetic  resonance  pictures  of  the  body  faster  and 
more  accurately). 

The  decomposition  of  signals  and  images  into  a  set  of  fundamental  building  blocks 
is  not  however  a  new  concept.  Fourier  analysis  is  such  a  decomposition.  Fourier 
methods  decompose  a  signal  into  a  series  of  sines  and  cosines.  The  discrete  set  of 
fundamental  frequencies  (the  “building  blocks”)  and  their  corresponding  amplitude 
completely  describe  the  initial  signal,  provided  some  care  is  taken  in  the  discrete 
sampling  procedure.  Various  properties  (e.g.  the  convolution  property)  combine  to 
make  Fourier  analysis  a  very  attractive  and  useful  mathematical  tool. 

Wavelets  are  however  an  entirely  new  set  of  building  blocks  that  present  a  number 
of  new  features.  One  advantage  of  the  new  wavelet  theory  lies  in  its  ability  to  zoom 
in  on  details  -much  like  a  camera  with  a  zoom  lens-  by  probing  signals  at  different 
scales  of  resolution.  Yet  another  advantage  lies  in  the  digital  processing  of  signals 
and  images:  by  processing  the  elementary  building  blocks,  instead  of  the  signal  or 
the  image,  the  processing  part  is  simplified,  and  new  insights  are  gained. 

A  novel  parametric  expansion,  called  Hermite-Rodriguez  expansion,  has  been  in 
troduced  as  functional  approximation  to  a  new  class  of  highly  singular  functions: 
hyperdistributions.  The  approximation  is  similar  to  the  way  a  Gaussian  of  narrowing 
width  is  a  Lighthill  approximation[lO]  to  the  Dirac  delta  “function’'  <5(x);  that  is,  the 
limit  of  the  sequence  of  Gaussians  of  narrowing  width  is  the  Dirac  delta  function: 

-xJ/A2 

lim  -  =.—  =  6(x)  (123) 

A—G 

We  will  now  prove  that  the  underlying  basis  functions  which  generate  Hermite- 
Rodriguez  expansions  are  in  fact  wavelet  functions[27). 

The  HR  expansion  of  a  function  f(x)  is  written  as: 

f(x)  =  £  anVnSx(x)  (124) 

n=0 

The  Rodriguez  formula  for  Hermite  polynomials[2]  is  then  used  to  obtain: 

OO 

/(x)  =  £an//„(*/A)<M*)  (125) 

n— 0 
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with 


■“F3C'lx,J,-(*/A>* 


(126) 


where  //n(ar)  denotes  the  Hermite  Polynomial  in  x  of  order  n.  We  now  redistribute 
2"  and  n!  in  the  following  way: 

/(x)  =  Q„ /W^(i  )  ^a(x)  where  an  =  /  Hdx(y)f{y)dy  (127) 

_ n  •/-O0 


where 


Hdl{x)  = 


Hn(x/\) 


(128) 


nV  '  2"/^ 

By  using  a  well-known[3]  bound  on  Hermite  Polynomials,  it  can  be  shown  that  the 
Hd  Polynomials  are  bounded  in  n: 

|  Hdx(x)  |<  Are*2/2*2  where  k  ~  1.086435  (129) 

Furthermore,  the  Hd  Polynomials  are  orthogonal  in  the  appropriate  L2  space,  with 


-*2/A2 


,  .  r+°°  .  ,  e~* 

<  Hdl  Hdxm  >  =  J  Hdtir)  «&(*)  dx  =  s/ii, 


(130) 


where  6nTn  denotes  Kronecker’s  delta;  the  generative  recursion  relation  for  Hd  Poly¬ 
nomials  becomes: 


Hdxn(x)  =  J~Hdxn_l(x)~ 


n  —  1 


with  Ild^(x)  =  1  and  IIdx(x)  =  — —  x 


(131) 


The  following  “analyzing” [27]  smooth  point- functions  are  now  introduced: 

Wl(*)  =  Hd*(z)6x{x)  (132) 

These  functions  present  the  following  properties: 

1.  Wx(  x )  is  a  function  of  class  C". 

2.  VF^(x)  has  a  rapid  Gaussian-like  decay  at  infinity.  In  other  words,  the  support 
of  the  function  is  said  to  be  essentially  compact. 

is  called  the  Hermite-Rodriguez  (HR)  wavelct[2l)  of  order  n  and  weight  X. 
The  Hermite- Rodriguez  (HR)  expansion  of  a  signal  f(x)  is  given  by 

00  r-f  00 

/(x)=5>nW*(x)  where  a„  =  /  Hdx(x)f(x)dx 

_ J -OO 


(133) 


A  sample  of  the  set  of  HR  wavelets  {W^x)},,  with  A  =  2  are  displayed  in  figure  5. 

In  a  similar  fashion,  the  HR  expansion  of  an  image  I(x,  y)  with  weights  A  and  [l 
is  given  by 

00  °°  r  r-f-00 

I(x,y)  =  J2  “nmW^fay)  where  anm  =  If  Hd*(x)Hd^{y)I(x,y)  dxdy 
n=0m=0  JJ-00 

(134) 

with  a  similar  convergence  condition  on  weights  (A,  jx),  and  with  the  2D  HR  wavelet 
function  defined  as  the  tensor  product  of  the  two  ID  independent  HR  wavelet  func¬ 
tions: 

=  W!(x)WZ(y)  (135) 

A  sample  of  the  positive  part  of  the  set  y)}n,n  for  A  =  2  is  displayed  in 

figure  6,  with  grey  scales  replacing  the  ordinate  values. 

Hermite- Rodriguez  wavelet  expansions  preserve  the  classical  properties  derivable 
for  orthogonal  polynomial  expansions  (as  opposed  to  power  series),  but  also  add  an 
important  new  feature:  the  weight  parameter  A.  The  convergence  condition  derived 
previously  from  Christoffel-Darboux  theory  holds  generally  on  a  semi-infinite  range 
of  A.  This  freedom  can  be  utilized  to  optimize  the  rate  of  convergence  of  the  HR 
wavelet  expansion.  HR  wavelet  expansions  are  also  found  to  be  robust  with  respect 
to  numerical  errors  on  the  HR  moment  coefficients.  In  other  words,  if  small  errors  are 
committed  on  the  exact  values  of  the  moment  coefficients  (for  example  by  quantizing 
the  moments  from  double-precision  reals  to  reals  with  one  or  two  decimal  places,  or 
even  to  integers),  the  convergence  properties  of  the  resulting  wavelet  expansion  do 
not  suffer. 


3.2  Hermite-Rodriguez  Wavelet  Expansion  of  a  Gaussian 

The  Hermite-Rodriguez  wavelet  expansion  of  a  Gaussian  of  unit  width  is  given  by: 

(136) 


n^O 

A  proof  based  on  linear  operator  relations  is  as  follows.  It  can  be  checked  by  differ¬ 
entiation  in  t  and  x  that 


=  etvi6(x) 


y/Ant 


Letting 
We  can  write 


At  =  A2 


-i2 /A2 

e  =  cx2v''46(x)  =  6x(x) 


v/jrA 


(137) 

(138) 

(139) 
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Consider  the  identity 


=  e^V’^M) 

=  e  ’  •’  y36A(x) 


Expanding  the  right-hand  side. 


evl/46(x) 


e-Vv^=E 

n=0 


(i  -  a 2r 

23r*/2\/nT 


»&(*) 


(140) 


(141) 


This  method  of  proof  provides  a  useful  check  on  the  numerical  calculations  of  HR 
moments. 


3.3  Hermite- Rodriguez  Wavelet  Expansions  in  One  Dimen¬ 
sion 


In  this  section,  the  convergence  properties  of  Hermite- Rodriguez  wavelet  expansions 
are  investigated.  To  this  end,  we  introduce  a  set  of  80  artificial  “landscapes”,  the  first 
16  of  which  are  displayed  in  figure  7.  These  landscapes  are  generated  by  summing  5 
Gaussians  whose  widths  and  offsets  are  generated  randomly  in  the  range  [0.1,1]  for 
the  widths,  and  [-3,3]  for  the  offsets.  The  landscapes  are  plotted  roughly  within  the 
support  of  three  standard  deviations,  ie.  in  the  range  [-6,6].  In  other  words,  if  Ln 
denotes  the  nth  landscape,  then: 

Ln(x)  =  J^6i'(x-xt)  (142) 

•=i 

where  <5,\(x)  denotes  the  Gaussian  of  width  A,  and  (A,,x,)  are  computer-generated 
pseudo-random  variables  with  respective  ranges  [0.1,1]  and  [-3,3].  All  landscapes  are 
then  rescaled  to  unit  height.  Each  landscape  is  sampled  by  501  datapoints  (501 
numbers  in  the  range  [-6,6]). 

If  all  5  widths  or  all  5  offsets  are  close  together  in  their  respective  values,  the 
associated  landscape  is  labeled  “realistic”  (ie.  the  landscape  resembles  a  realistic 
mountain  or  mountain  range),  while  if  on  the  contrary,  the  5  widths  and  offsets  are 
spread  out  in  range,  the  generated  landscape  is  labeled  “futuristic”. 


Although  Hermite  polynomials  of  high  order  are  known  to  diverge  polynomially 
(xn,  “whipping-tail”)  as  their  order  n  increases,  the  Hermite- Rodriguez  wavelets  of 
high  order  are  well-behaved,  and  decreasing  in  modulus.  More  specifically,  a  classical 
bound  of  Hermite  Polynomials[3]  yields: 


I  KM  l<  k 


€~Z^2XJ 

x/jrA 


where  k  m  1.086436 


(143) 
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Figure  8  displays  the  HR  wavelets  of  orders  1000  and  10000. 

In  Figure  9,  a  shifted  Gaussian  (centered  at  x  =  2),  of  width  A  =  1/2,  renormalized 
to  unit  height,  denoted  by  Lo(x ),  is  approximated  with  a  Hermite-Rodriguez  wavelet 
expansion  A\(x)  of  weight  A  =  1,  with  Ntk  partial  sum: 

N  ,+6 

^4^(z)  =  ^2  cii  Wflx)  where  a,  =  /  L0(x)Hdf(x)  dx  (144) 

i=0 

where 

W?(*)  = /«?(*)«»(*)  and  //<(?( (145) 

2n/2Vn! 

and  Hi(x)  denotes  the  Hermite  polynomial  of  order  i.  The  expansion  is  seen  to 
converge.  The  40*A  partial  sum  already  closely  approximates  the  original  Gaussian. 

In  figure  10,  a  futuristic- looking  landscape  (landscape  #68  and  denoted  by  Les{x)) 
is  approximated  with  a  Hermite-Rodriguez  wavelet  expansion  of  width  A  =  1.  The 
300*/l  partial  sum  closely  approximates  the  original  landscape.  A  high  number  of 
partial  sums  is  required  for  the  capture  of  high  frequency  contents. 

In  figures  11  and  12,  a  less  futuristic-looking  (but  still  futuristic)  landscape  (land¬ 
scape  #71  and  denoted  by  L7i(x))  is  approximated  with  a  Hermite-Rodriguez  wavelet 
expansion  for  two  different  values  of  the  weight  A:  1.0  and  1.15  respectively.  The 
140'/l  partial  sum  for  A  =  1.15  quickly  approximates  the  original  landscape,  while 
the  wavelet  expansion  for  A  =  1.0  has  a  much  slower  convergence  rate.  This  free 
parameter,  can  thus  be  utilized  to  improve  rate  of  convergence.  Figures  13  through 
17  emphasize  the  critical  role  of  the  weight  parameter  A  in  improving  the  conver¬ 
gence  of  Hermite-Rodriguez  wavelet  expansions.  Figure  18  displays  the  robustness  of 
HR  expansions,  when  HR  moments  are  rounded  off.  The  rate  of  convergence  of  the 
quantized  wavelet  expansion  is  only  very  slightly  modified. 

It  can  be  observed  that  landscapes  with  an  increasing  number  of  higher  frequencies 
require  an  increasing  number  of  partial  sums  for  accurate  HR  wavelet  expansion 
representations.  This  was  expected  for  HR  wavelet  expansions,  which  are  of  global 
nature  (moment  expansions  vs.  Taylor  expansions),  and  capture  global  behavior  very 
fast,  while  they  are  slower  in  capturing  sharply  localized  features. 
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3.4  Hermite- Rodriguez  Wavelet  Expansions  in  Two  Dimen¬ 
sions 


In  a  similar  fashion,  the  HR  wavelet  expansion  of  an  image  I{x,y)  with  weights  A 
and  p  is  given  by 

00  00  ,  rt+oo 

=  E  E  where  «nm  =  //  Hd*(x) Hd^{y) I (x , y)dxdy 

n=0m=0  JJ~°° 

(H6) 

with  the  2D  HR  Wavelet  function  defined  as  the  tensor  product  of  the  two  ID  inde¬ 
pendent  HR  Wavelet  functions: 

wlt‘‘\z.y)  =  w^x)WMy)  (147) 

A  sample  of  the  positive  part  of  the  set  for  A  =  p  =  2  is  displayed 

in  figure  6,  with  grey  scales  replacing  the  ordinate  values. 


A  set  of  SO  artificial  “nebulae”,  the  first  16  of  which  are  displayed  in  figure  19. 
These  nebulae  are  generated  by  summing  5  2D  Gaussians  whose  widths  and  offsets 
are  generated  randomly  in  the  range  [0.2,1]  x  [0.2,1]  for  the  widths,  and  [-2,2]  x  [-2,2] 
for  the  offsets.  The  nebulae  are  plotted  within  the  range  [-6,6]  x  [-6,6].  In  other  words, 
if  Nn  denotes  the  nth  nebula,  we  have: 

5 

Nn(x)  =  ~  xt,y  -  y>)  (MS) 

•=i 


where  (x,  y)  denotes  the  2D  Gaussian  of  width  A;  and  /*,,  which  is  the  tensor 
product  of  two  ID  Gaussians: 


(M9) 


and  (A,,//,)  and  (x,,y,)  are  computer-generated  pseudo-random  variables  with  respec¬ 
tive  ranges  [0.2,1]  and  [-2,2].  All  nebulae  are  then  rescaled  to  unit  height.  Each  nebula 
is  sampled  by  64x64  datapoints. 


4  Green’s  Function  for  the  AntidifFusion  Equa¬ 
tion 

In  this  section,  we  repeat  in  more  detail  two  proofs  given  in  sections  1  and  2.  Namely 
that  Fourier  transforms  do  not  provide  us  with  a  Green’s  function  for  the  antidiffusion 
equation  (i.e.  a  convolution  inverse  of  a  gaussian),  but  that,  in  contrast,  there  is  a 
simple  hyperdistribution  representation  for  that  same  Green’s  function.  We  then 
attempt  to  reconstruct  images  corrupted  with  gaussian  blur  by  letting  the  blurred 
images  evolve  freely  as  initial  conditions  of  the  antidiffusion  equation. 
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4.1  Antidiffusion  and  Fourier  Transforms 


Consider  the  Gaussian  filter  G(x,t)\ 


G{x,t)  = 


v/4  nt 


(150) 


We  suppress  the  spatial  variable  x  if  confusion  does  not  arise,  and  denote  G(x,t) 
by  G(t).  We  have  shown  in  section  1.2  that  an  attempt  to  evaluate  an  explicit 
representation  of  the  convolution  inverse  of  G(t)  from  Fourier  transforms  yields: 


1  f+a 

F(a,T,2t)  =  -/_ 


-ikx+k2t 


i^[erf(  iVi  a+irt)+ Crf(,v' 0 "  in)] 


(151) 


(152) 


where  erf(x)  denotes  the  error  function[3].  Proceeding  formally, 

7nu[G(f)]  =  lima_00  F{a,x,2t),  which  was  shown  to  diverge. 

G(t)  is  a  smooth  point  function  for  nonegative  values  of  t.  By  contrast,  Inv{G{t)\  is 
not  even  a  distribution,  since  it  does  not  exhibit  weak  convergence^]- [1 1].  Specifically, 

i2  . 

its  integral  over  the  Gaussian  test  function  e”  «  diverges.  Proceeding  with  the  use  of 
Fubini’s  theorem[25]  on  the  interchange  of  the  limit  and  integral  sign, 

/+°o  r+a  J  r+oo  _  *2  , — 

e  i««  lim  F(a,x,2<)  dx  =  lim  /  e  c  41  '  1  dx  dk  =  lim4aV7rf  — *  +oo 

-oo  a— 00  a— oo  J_a  J_(X>  »— oo 

(153) 

In  other  words,  the  convolution  inverse  of  a  Gaussian  filter  as  given  by  Fourier  trans¬ 
form  methods  is  not  a  smooth  point-function. 


4.2  Antidiffusion  and  Hyperdistributions 

In  the  previous  paragraph,  it  was  proven  that  Inv[G{t)\,  as  obtained  by  Fourier 
transform  methods,  is  ill-defined.  By  contrast,  we  can  give  an  explicit  representa 
tion  of  Inv\G{t)}  using  Lyperdistributions.  The  Green’s  function  for  the  diffusion 
equation[21][22]  can  be  written[23]  as 

G(x,<)  =  etv2S(x)  (154) 

where  6  denotes  the  Dirac  delta  function,  and  which,  for  t  >  0,  converges  to  the 
well  known  Gaussian  given  above  in  Eq.(150).  The  function  G ,  also  called  the  initial 
condition  Green’s  /uncfion[24],  has  the  properties 

^G(x,t)=V2G(xJ)  ,  G(x, 0)  =  fi(x)  (1 55) 

at 
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We  note  that  the  two  infinite  series  of  distributions: 


e^’six)  =  -£  i(±  (156) 

are  hyperdistributions.  Furthermore,  the  expressions  are  inverse  to  each  other  under 
convolution,  as  can  be  checked  by  the  Bochner  algebra  for  hyperdistributions.  The 
convolution  inverse  Inv[G(t)\  of  the  Gaussian  is  now  readily  obtained  using  Eq.(154). 
Thus 

6(x)  =  (e‘v^(x))  *  7no[G(t)]  =  e^\S  *  /nu[G(t)])(x)  =  etV*  Inv[G{t)]  (157) 
from  which  we  conclude  that 

/nu^O]  =  e-‘y2<5(x)  (158) 

We  note  the  relation  of  the  convolution  inverse  Inv[G(t)\  to  the  time  reflected  solution 
of  the  diffusion  equation: 


Inv[G{t)\  =  G(-t)  =  GT(t )  (159) 

where  T  denotes  ”time  reflected”.  The  quantity  G T(t)  satisfies  the  equation: 

jtaTw  =  -vjgt(<)  ,  gt(o)  =  %)  C60) 

and  therefore  /nt>[G(f)]  is  the  fundamental  solution  of  the  antidiffusion  equation 
(time-reversed  diffusion).  The  sequence  {F(a,x,2t)}„  plotted  in  figure  4  can  be  shown 
to  define  the  same  Hyperdistribution  given  by  Eqs.(156)  and  (158).  The  sequence  in 
Eq.(156)  is  analogous  to  the  sequence  that  defines  8(x)  in  terms  of  narrowing  Gaus- 
sians.  Both  G  and  Inv[G ]  have  analytic  power  series  in  wave-number  space  with 
Eq.(156)  as  Fourier  images.  This  is  not  accidental;  in  fact,  hyperdistributions  have 
the  algebraic  properties  of  formal  power  series,  since  the  Fourier  transform  of  hyper¬ 
distributions  is  a  power  series  in  the  wave-number  k. 


4.3  Antidiffusing  Images  Corrupted  with  Gaussian  Blur 

We  have  shown  in  the  previous  section  that  the  Bochner  algebra  for  deconvolving  hy¬ 
perdistributions  leads  to  a  theoretically  possible  reconstruction  of  waveforms  blurred 
by  a  Gaussian  filter  by  antidiffusing  the  blurred  waveform  for  an  appropriate  duration 
of  time.  In  this  section,  we  will  go  through  the  proof  in  more  detail. 

Suppose  an  image  is  blurred  by  convolution  with  a  two-dimensional  Gaussian  filter 
of  width  A.  The  convolution  inverse  to  the  Gaussian  filter  is  obtained  with  the  help 
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of  the  Bochner  algebra  for  hyperdistributions.  The  2D  Hermite-Rodriguez  moments 
of  a  normalized  Gaussian  filter  of  width  A  are  given  exactly  by 


ff+oo  e  —Jr 

anm  =  JJ-  oo  v^A 


2  2 

— x  -iJL  (  n'm*  *r  j 

xn  ym  dxdy=\  (f )!(?)!  *3<-+")/*  lfnandrr 
V’rA  I  0  otherwise. 


m  are  even. 


Thus,  the  hyperdistribution  expansion  of  the  Gaussian  filter: 

+  3Q  +oo  / _ 1  / _ I  \m 

h(x,y)  =  E  E 


n=0  m=0 


+oo  +oo  \4  j 

=  E  E< j)"+"^rvrv;-<(x)«(!,) 

n=0m=0  °  Uin. 


which  can  be  separated  in  the  form 


The  above  double  m  can  be  written  in  simple  operator  form  as 
6\(x,y)  =  e~i~v*6(x)e~»'v*6(y)  =  e"*“v2^(ar)^(j/) 


(161) 


(162) 


«*■»>  =  {Sir^*"^)}  {|0(f  )-sr^(»)}  <163> 


(164) 


The  convolution  inverse  /nu[<5\]  of  <5.*  is  now  easily  determined  through  the  2D 
Bochner  Algebra.  More  specifically,  it  represents  a  special  simplifying  case  where  the 
2D  algebra  reduces  to  the  product  of  two  ID  Bochner  algebras.  That  is  because  the 
above  hyperdistribution  of  two  independent  variables  breaks  down  into  the  product 
of  two  hyperdistributions  of  a  single  variable  (respectively  x  and  y)  as  can  be  seen  in 
Eq.(163).  This  is  not  surprising;  a  two-dimensional  Gaussian  is  separable  in  its  two 
independent  variables.  Thus,  for  all  sujh  “radial”  filters,  the  Bochner  algebra  reduces 
to  the  product  of  two  ID  algebras: 


Cnm  =  |l>n-p&pj  j£an-P6pj  (165) 

where  the  superscripts  1  and  2  refer  to  the  respective  independent  variables  x  and  y. 
The  convolution  inverse  to  6\  is  then  easily  determined  (see  appendix  A)  to  be 


+oo  +oo  \  4  1 

_ n  _n  O  Tl.TTl. 

n=u  m=u 


(166) 


which  can  be  written  again  as 


/nu[6A](x,i/)  =  e~  *•  x*6{x,y) 


(167) 


31 


Thus,  filtering  an  image  7(x,y)  with  6*(x,y)  yields  the  blurred  image  J(x,y ) 
where 

J(x,y)  =  8\(x,y)  *  I(x,y)  =  eVv\5(x,y)  *  7(x,y)  =  eVyJ/(i,y)  (168) 
And  finally,  reconstructing  the  image  /{x,  y): 

7(x,y)  =  J(x,y)  *  /nn[5>](i,y)  =  e"VyJ£(x>y)  *  J(i,y)  =  e"VyJ  J(x,y)  (169) 

The  use  of  linear  operators  allows  simple  representations  for  all  the  quantities  in¬ 
volved. 

Eq.(169)  can  then  be  further  simplified  with  the  help  of  Taylor  expansions  in  the 
following  way.  Write  down  A 4/8  as  a  finite  sum  of  equally  small  quantities.  Then 

A4/S  =  E*«o  and 

7(x,y)  =  e-Vv2J(x,y)  =  J(x,y)  =  J]  {e"A*v2 J(x,y)}  (170) 

And  since  A*  is  small,  e  -A*yJ  can  be  approximated  to  first  order  by  its  Taylor  ex 
pansion: 

e~A‘vJ  =  1  -  A*V2  +  0{&l)  (171) 

Eq.(170)  can  thus  be  approximated  to  first  order  by: 

'(*.y)  =  n  {•/(*,!,)- A  VV(*,j,)}  (172) 

where  the  subscript  k  in  A  has  been  dropped.  It  can  now  be  shown  that  the  origi¬ 
nal  image  can  be  recovered  by  running  a  iime-reversed  diffusion  with  the  degraded 
image  J  as  initial  condition.  The  normalized  time-reversed  diffusion  equation  (or 
antidiffusion)  can  be  written  as 

-J  +  V3J  =  0  (173) 

ot 

and  can  be  forward-discretized  in  time  as 

=  -V2J(  (174) 

or  equivalently  as: 

Jt+*t  =  Jt  ~  AW,  (175) 

Equations  (172)  and  (175)  are  one  and  the  same,  and  therefore  the  original  image  can 
theoretically  be  reconstructed  by  applying  a  time-reversed  diffusion  on  the  degraded 
image  with  a  very  small  time-step  (so  that  Eq. ( 171)  holds  to  first  order).  Higher-order 
corrections  to  our  approximation  can  be  implemented  by  expanding  the  exponential 
in  Eq.(170)  to  higher  order. 
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We  dow  discretize  the  2+1  antidiffusion/diffusion  equation  with  diffusivity  ±v: 


±  i/V2u  =  0 


(176) 


in  the  following  way.  The  infinite  plane  is  sampled  at  locations  on  a  square  grid: 


U  - >  U: 


(177) 


where  the  subscripts  i  and  j  represent  the  spatial  location  of  a  mesh  point,  and  the 
superscript  k  represents  the  time  index.  The  continuous  operators  V2  and  d/dt  are 
replaced  by  the  discrete  approximations: 


(178) 


where  St  represents  the  discretized  time-step,  (this  is  a  forward-in-time  discretization 
scheme,  with  an  error  of  the  order  of  S 2),  and 


2  k  QlVlo  +/?Eu*|o  +7«fj 

(q  +  0)6} 

where  6X  represents  the  discretized  space-gap,  and 


(179) 


7  =  -4(0  +  /?)  (180) 

uk  |o  represents  the  sum  of  all  immediate  neighbors  of  u*  forming  a  diamond  pat¬ 
tern  centered  at  it,*,  and  J2U>C  |o  represents  the  sum  of  all  immediate  neighbors  of 
ukj  forming  a  square  pattern  centered  at  uk-  (this  is  a  central-in-space  discretization 
scheme,  with  an  error  of  the  order  of  <^).  The  configuration  of  neighbors  forming  a 
diamond  pattern  is  called  the  diamond  stencil ,  while  the  configuration  of  neighbors 
forming  a  square  pattern  is  called  the  square  stencil.  This  central-in-space  discretiza¬ 
tion  of  the  continuous  operator  V2  is  carried  out  in  appendix  B.  How  is  the  ratio 
a//3  to  be  chosen[28]?  In  other  words,  what  is  the  weighting  of  the  diagonal  and  the 
square  neighbors  in  the  discretization  of  the  operator  V2?  It  is  found  that  the  best 
scheme  involves  weighing  diamond  and  square  stencils  based  on  their  mesh-distance 
from  the  central  grid-point  (square  neighbors  are  by  a  factor  y/2  further  apart  from 
the  central  grid- point  than  diamond  neighbors): 

a//?  =  \/2  (181) 


The  numerical  values  chosen  in  the  discretization  are  the  following.  6t  must  be 
small  for  the  approximation  in  Eq(171)  to  hold  to  second  order,  and  ST  must  be  small 
for  the  approximation  in  Eq.(179)  to  hold  to  fourth  order.  More  specifically,  we  have 
chosen  v  =  1  L2/T,  8X  =  12/65  A,  8t  =  1.7  x  10~6T,  where  L  is  a  unit  of  length. 
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and  T  is  a  unit  of  time.  The  values  of  t  specified  in  the  figures  have  been  rescaled  to 
higher  values  for  ease  of  comparison.  In  order  to  relate  these  diffusion  durations  to 
the  duration  of  time  T  responsible  for  the  diffusion  eTV*  I(x,y)  of  the  image  I(x,y), 
divide  t  by  the  factor  0.0085. 

In  figures  20  through  23,  images  of  the  letters  “T”  and  “E”  are  diffused,  and 
a  reconstruction  is  attempted  by  antidiffusing  the  blurred  images.  If  the  diffusion 
duration  is  too  long,  the  attempted  reconstruction  fails. 

In  figure  24,  an  initial  pattern  in  the  shape  of  a  spiral  is  diffused,  and  a  recon¬ 
struction  is  attempted  for  increasing  amounts  of  diffusion  duration.  Once  again,  it  is 
observed  that  if  the  diffusion  duration  is  too  long,  the  reconstruction  fails. 

In  figure  25,  the  effect  of  noise  on  reconstruction  by  antidiffusion  is  analyzed.  An 
intial  pattern  in  the  shape  of  a  letter  “T”  is  diffused  for  a  certain  amount  of  time.  A 
2D  field  of  pseudo-white  noise  is  generated  with  a  random-number-generator.  That 
field  is  subsequently  added  onto  the  diffused  pattern  at  a  proportion  given  by  the 
Signal-to-Noise  Ratio  (SNR).  The  resulting  pattern  is  subsequently  antidiffused  for 
the  appropriate  amount  of  time,  and  the  resulting  reconstructions  are  displayed. 

We  conclude  by  noticing  that  the  reconstruction  by  antidiffusion  of  blurred  images 
is  an  unstable  algorithm,  with  the  noise-content  of  the  image  (or  numerical  round-off 
errors)  as  the  source  of  the  instability.  Each  mode  is  exponentially  amplified  in  mag¬ 
nitude  with  antidiffusion.  The  Signal-to-Noise  Ratio  thus  plays  a  critical  perturbative 
role.  The  reconstruction  process  breaks  down  if  the  width  of  the  Gaussian  filter  is 
too  wide  (i.e.  if  the  diffusion  and  thus  antidiffusion  durations  are  too  long:  figures 
21,  23,  and  24),  or  if  the  Signal-to-Noise  Ratio  is  too  low5  (figure  25). 

In  the  next  section,  we  will  approximate  hyperdistribution  sums  with  Hermite- 
Rodriguez  Wavelets,  and  once  again  attempt  to  reconstruct  signals  and  images  cor¬ 
rupted  with  gaussian  blur,  in  the  difficult  case  of  low  Signal-to-Noise  Ratios.  We  will 
then  generalize  my  results  to  include  the  case  of  non-gaussian  blur. 


5  Deblurring  with  Hermite-Rodriguez  Wavelets 

5.1  Introduction 

In  the  previous  section,  we  have  concluded  that  deblurring  by  discretization  of  the 
antidiffusion  equation  is  a  chaotic  process,  since  numerical  errors  increase  exponen¬ 
tially  with  time.  In  other  words,  antidiffusion  is  extremely  sensitive  to  noise,  and  an 

5  round-off  error  in  the  computation  is  equivalent  to  noise  added  onto  the  blurred  image.  The  com¬ 
putation  is  implemented  with  double-precision  fortran  real  variables,  corresponding  to  a  phenomenal 
1015  SNR,  as  long  as  the  blurred  image  is  free  of  any  other  kind  of  noise 
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extremely  high  Signal-to-Noise  Ratio  (SNR)  is  necessary  for  reconstructing  a  blurred 
image  via  antidiffusion[14].  In  this  section,  we  modify  the  antidiffusion  algorithm  by 
expanding  the  blurred  waveform  in  the  form  of  a  Hermite- Rodriguez  wavelet  expan¬ 
sion.  The  antidiffusion  is  then  carried  out  on  the  building  blocks  of  the  expansion  -the 
HR  wavelets-,  as  opposed  to  the  waveform  as  a  whole.  Hermite-Rodriguez  wavelets 
are  known  exactly,  and  can  thus  be  more  successfully  antidiffused  when  the  blurred 
waveform  is  corrupted  by  noise.  The  noise  enters  effectively  only  on  the  Hermite- 
Rodriguez  moments,  which  are  constant  and  remain  unchanged  by  antidiffusion.  This 
deblurring  procedure  is  found  to  be  more  robust  with  respect  to  noise.  Artificial  ID 
signals  (“landscapes”)  and  artificial  2D  astronomical  images  (“nebulae”),  which  have 
been  blurred  by  a  diffusion  process  and  subsequently  corrupted  with  white  noise,  are 
successfully  reconstructed. 

5.2  Hermite-Rodriguez  Wavelets  and  Antidiffusion 

The  fundamental  solution  of  the  diffusion  equation [21] [22]  at  time  t  can  be  written[23] 
as 

G(x,f)  =  e‘v26(x)  (182) 

where  S  denotes  the  Dirac  delta  function,  and  which,  for  t  >  0,  converges  to  the 
Gaussian  form: 

g(i,<)=7aT  (183) 

The  function  G,  also  called  the  initial  condition  Green’s  function[2A],  has  the  proper¬ 
ties 

§-G(x,t)  =  V2G(x,f)  ,  G(x,0)  =  «5(x)  (184) 

It  has  been  proven  that  the  two  hyperdistributions  defined  by 

e±re’i(i)  =  f;  i(±  1  V!)"<S( x)  (185) 

n=:0 

are  inverse  to  each  other  under  convolution.  That  is,  the  convolution  e~tVl 6{x)  * 
e+ty2£(x)  yields  the  Dirac  delta  function  S(x).  Consequently,  the  fundamental  solu¬ 
tion  of  the  antidiffusion  equation 

^G(x,0  +  V2C(x,0  =  0  (186) 

at 

can  formally  be  written  as  e~*y2£(x). 

The  set  of  HR  Wavelets  |  W^(x)}a  can  be  shown  to  be  a  globally  invariant  subspace 
(in  A)  of  the  operators  e±<v  (see  appendix  C).  In  other  words,  it  can  be  proven  that, 
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diffusing  the  HR  wavelet  of  order  n  and  weight  Aj  for  a  duration  of  time  t  yields  a 
rescaled  HR  wavelet  of  same  order  but  higher  weight  A2: 


,,+tv2 


/Vf 

[hj 


(187) 


where  A2  =  y  Aj  +  4 1.  Consequently,  the  above  statement  can  be  time-reversed,  to 
show  that  the  antidiffusion  of  the  HR  wavelet  of  order  n  and  weight  A2  for  a  duration 
of  time  t  yields  a  rescaled  HR  wavelet  of  same  order  but  lesser  weight  Aj: 

Wt(x)  (188) 

This  property  is  utilized  in  the  following  section. 


5.3  Wavelet  Reconstruction  of  Diffused  Signals  and  Images 
(Gaussian  Blur) 


A  one-dimensional  signal  S(x )  is  diffused  for  a  duration  of  time  t.  The  resulting 
blurred6  signal  is  denoted  as  S(x).  The  goal  is  to  reconstruct  S(x)  by  antidiffusing 
5(x)  for  the  same  duration  of  time  i.  To  this  end,  we  expand  5(x)  as  the  Nth  partial 
sum  of  an  HR  expansion  of  weight  A2: 

W*(x)  (189) 

n=0 


where  {d„}„  are  the  Hermite  moments  of  the  blurred  signal  S{x),  weighted  by  A2. 
The  partial  sum  above  can  be  formally  antidiffused  by  antidiffusing  the  constitutive 
blocks  of  the  expansion,  which  are  the  HR  wavelets:  by  formally  antidifTusing  for  a 
duration  of  time  t,  we  obtain: 


e~tV7S{x) 


,-tv2 


(Ea-vet*))  =  =  E*.  (¥] 

n=0  n=0  n=0  \A1/ 

(190) 

with  A]  =  \J A|  —  4t.  The  antidiffused  signal  can  therefore  be  put  in  a  form  very 
similar  to  an  HR  expansion,  this  time  of  weight  Aj.  However,  the  antidiffused  signal 
in  Eq.(190)  is  not  an  HR  expansion,  because  of  the  extra  factor  multiplying  the 

Hermite  moments  an.  Consequently,  a  condition  of  the  same  type  as  in  Eq.(l  14)  does 
not  guarantee  the  convergence  of  the  antidiffused  series.  This  matter  will  be  treated 
in  the  next  paragraph.  The  case  of  two-dimensional  images  is  similar,  since  two- 
dimensional  HR  wavelets  are  the  tensor  product  of  two  one-dimensional  HR  wavelets. 


5in  this  section,  we  restrict  ourselves  to  pure  gaussian  blur 
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As  a  result,  the  antidiffused  HR  expansion  of  an  image  can  be  put  in  a  form  similar  to 

(,  \  n+m 

j*J  multiplying  the  Hermite  moment 

of  order  n  x  m. 

Blurred  signals  and  images  (pure  gaussian  blur)  have  been  formally  reconstructed 
by  antidiffusing  their  HR  expansions.  Since  e±tVl  is  a  linear  operator,  it  is  known[29] 
that7  if  both  series 


S(x)  = 
R(x)  = 


(191) 

(192) 


are  convergent,  then  the  reconstructed  signal  is  given  by  R(x )  =  e~tV* S(x),  with 
R(x)  given  by  Eq.(192).  However,  there  is  no  guarantee  that  the  series  R(x)  does 
converge. 


More  specifically,  if  the  moments  d„  do  not  decrease  to  zero  as  fast  as 
diverges  a s  n  — ►  oo,  then  the  series  R(x)  diverges.  This  is  in  fact  the  typical  case 
when  the  diffused  waveforms  are  corrupted  by  instrumental,  numerical,  or  computer 
round-off  noise.  To  see  how  this  happens,  we  decompose  the  Hermite  moments  an 
for  a  diffused  and  noisy  signal  as  the  sum  of  the  analytically  exact  (but  unknown) 
moments  bn  ,  and  the  residual  error  en: 


«n  —  bn  +  (193) 

The  antidiffusion  of  the  HR  expansion  with  the  analytically  exact  moments  bn  yields 
precisely  the  HR  expansion  which  converges  to  the  original  signal  S(x): 

{oo  00 

£  *»»?>(*)  i  =  £  U  r  lC(*)  =  £  a„W„A'(x)  (194) 

n=0  )  n=0  \/'1/  n=0 

On  the  other  hand,  the  antidiffusion  of  the  HR  expansion  with  the  residual  error 
coefficients  en  diverges: 


({^)  W^'(i)  =  oo  (195) 

However,  since  usually  en  <C  6„,  the  series  in  Eq.(195)  will  require  a  minimum  amount 
of  partial  sums  before  the  residual  error  can  grow  in  modulus  and  overtake  the  exact 
solution  as  given  by  Eq.(194).  The  series  in  Eq.(195)  is  called  partially  convergent,  or 

‘in  one  dimension,  without  loss  of  generality 
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asymptotic.  Consequently,  if  the  Nth  partial  sum  EiL,  «.»•*>(*)  is  a  good  approxi¬ 
mation  of  the  original  signal  S(x ),  and 


(196) 

n=0  \  /  n=0 


which  usually  implies  that 


then  the  reconstruction 


(197) 


(198) 


is  successful,  and  R(x)  «  S(x).  If,  on  the  other  hand,  the  antidiffused  residual  error 
series  grows  too  fast  (e.g.  if  b„  and  e„  are  of  the  same  order,  which  is  the  case  if  the 
Signal-to-Noise  Ratio  is  too  low),  and  higher-order  partial  sums  are  required  for  good 
approximations  of  the  signals  5(i)  and  S( x)  (e.g.  if  the  original  signal  comprises 
a  large  number  of  high-frequency  components),  then  the  reconstruction  will  fail.  In 
other  words,  for  lower  order  partial  sums,  antidiffusion  can  be  implemented  for  longer 
durations  of  time,  and  consequently  reconstruct  severly  blurred  signals.  If,  on  the 
other  hand,  good  approximations  of  the  signals  5(ar)  and  S(x )  require  partial  sums 
of  higher  order,  the  amount  of  antidiffusion  which  can  be  performed  on  the  blurred 
signal  before  instrumental  and  numerical  errors  overcome  the  calculation  is  limited. 
This  limitation,  as  well  as  successful  deconvolutions,  are  displayed  in  figure  26  for 
one-dimensional  signals,  and  in  figure  27  for  two-dimensional  images. 


We  now  investigate  whether  it  is  possible  to  determine  a  priori  the  threshold  in  the 
order  N  for  which  the  reconstruction  by  antidiffusion  will  fail.  To  this  effect,  We  can 
utilize  one  of  two  classical  constraints  employed  in  signal  and  image  reconstruction: 


•  positivity:  the  intensity  of  each  pixel  of  the  reconstructed  waveform  must  be 
positive,  since  the  original  waveform  is  everywhere  positive. 

•  total  intensity  conservation:  the  total  intensity  of  the  waveform  is  a  con¬ 
served  quantity  under  any  diffusion/antidiffusion  process  (see  appendix  E).  Con¬ 
sequently,  it  is  required  of  the  reconstructed  waveform  to  conserve  the  total 
intensity  of  the  blurred  waveform. 


At  a  certain  threshold  in  the  order  N  of  the  reconstruction,  one  of  the  two  criteria 
above  is  destined  to  fail.  At  that  point,  the  reconstruction  will  fail  as  well.  That 
is  because  the  diverging  antidiffused  residual  error  series  in  Eq.(195)  overtakes  in 
magnitude  the  value  of  the  exact  antidiffused  series  in  Eq. (194).  The  upper  limit  in 
1 V  for  which  the  two  <~riteria  above  are  still  verified  will  yield  the  order  of  the  best 
possible  reconstruction  by  antidiffusion  of  the  HR  expansion  of  the  blurred  waveform. 
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5.4  Application  to  Simulated  Landscapes  and  Astronomical 
Images  Corrupted  with  Gaussian  Blur 

In  figures  26  and  27,  a  one-dimensional  signal  (landscape  11)  and  a  two-dimensional 
image  (nebula  62)  are  diffused  for  a  certain  amount  of  time.  The  diffused  waveforms 
are  subsequently  corrupted  with  white  noise  obtained  by  a  pseudo-random  number 
generating  algorithm.  The  one-dimensional  signal  is  corrupted  with  3%  additive  noise 
(i.e.,  the  noise  is  added  uniformly  onto  the  diffused  signal,  at  a  proportion  of  3%  of  the 
maximum  intensity  of  the  signal).  The  two-dimensional  image  is  corrupted  with  30% 
multiplicative  noise  (i.e.  the  noise  is  added  uniformly  onto  the  diffused  image,  pixel  by 
pixel,  at  a  varying  proportion  of  30%  of  the  intensity  of  each  pixel).  Both  the  diffused 
and  noisy  signal  and  image  are  then  approximated  by  a  ID  and  2D  HR  expansion 
respectively.  The  HR  expansions  are  subsequently  antidiffused  by  antidiffusing  the 
underlying  HR  wavelets.  The  antidiffusion  of  the  ID  HR  wavelet  of  order  n  and 
weight  A  for  an  amount  of  time  t  is  known  exactly:  it  yields  the  HR  wavelet  of  same 
order,  of  weight  y/X1  —  At,  rescaled  by  tue  factor  {A/(\/A2  —  4f)}n.  The  antidiffusion 
is  carried  out  on  the  blurred  signal  (for  the  same  amount  of  time  as  the  earlier  diffusion 
duration)  on  two  HR  partial  sums  of  different  order,  one  of  lower  order,  the  other 
of  higher  order.  Both  reconstructions  are  displayed.  The  antidiffusion  of  the  lower- 
order  partial  sum  is  successful,  while  the  antidiffusion  of  the  higher-order  partial  sum 
displays  the  diverging  effects  of  noise  on  the  Hermite  moments.  The  antidiffusion  is 
also  carried  out  on  the  blurred  image  for  two  different  HR  partial  sums.  The  same 
remarks  apply.  The  tests  of  positivity  and  conservation  of  total  intensity  are  carried 
out  for  the  signal  in  figure  26  and  the  image  in  figure  27.  The  results  are  displayed 
in  figure  28,  together  with  a  least-squares  goodness-of-fit  evaluation  of  the  fidelity  of 
the  reconstructed  waveform  with  the  original  waveform,  plotted  against  the  order  N 
of  the  partial  sum.  The  best  reconstruction  is  located  at  the  minimum  of  the  least- 
squares  curve,  and  this  does  correspond  to  the  upper  limit  in  the  order  N  for  which 
the  reconstructed  waveform  is  everywhere  positive,  and  for  which  the  total  intensity 
remains  equal  to  the  total  intensity  of  the  blurred  waveform. 

It  is  noteworthy  to  point  out  that  the  noisy  signals  and  images  were  not  pre- 
processed  for  noise  reduction  before  antidiffusion. 

In  figure  29,  the  blurring  of  the  image  in  figure  27  is  repeated,  and  a  reconstruc¬ 
tion  of  the  original  image  by  a  direct  discretization  of  the  antidiffusion  equation  is 
attempted.  In  the  case  when  the  blurring  is  clean  -that  is,  when  the  only  source  of 
noise  on  the  diffused  image  is  numerical  round-off  noise8-,  the  reconstruction  is  suc¬ 
cessful.  In  the  case  when  the  diffused  image  is  corrupted  by  noise,  the  reconstruction 
is  highly  unsuccessful. 

8the  discretization  is  performed  with  double-precision  real  variables,  bringing  the  Signal-to-Noise 
Ratio  in  this  case  to  about  1015. 
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The  advantages  of  wavelet-based  antidiffusion  as  an  algorithm  for  deblurring  sig¬ 
nals  and  images  are: 

•  The  anlidiffusion  algorithm  does  not  require  any  lengthy  computation,  since 
antidiffused  wavelets  are  known  exactly.  This  is  in  sharp  contrast  with  the 
antidiffusion  algorithm  (directly  involving  a  discretization  of  the  antidiffusion 
equation,  which  has  a  very  high  computational  load),  or  even  Fourier  transform 
methods  (which  require  lengthy  integral  calculations). 

•  The  algorithm  is  noise-robust,  since  noise  enters  only  in  the  computation  of  the 
Hermite  moments  of  the  blurred  waveform.  In  other  words,  noise  is  involved  in 
an  integration  (the  moment  computation),  which  is  a  “smoothing”  operation, 
as  opposed  to  a  differentiation  process  (e.g.  direct  discretization  of  the  antid¬ 
iffusion  equation),  which  is  a  “roughening”  operation.  Nevertheless,  noise  or 
numerical  errors  on  the  Hermite  moments  eventually  impair  the  reconstruction 
of  the  blurred  waveform  in  the  case  where  higher-order  partial  sums  are  nec¬ 
essary  for  accurate  HR-expansion  approximations  of  the  blurred  and  original 
waveforms,  as  well  as  in  the  case  of  large  diffusion  durations. 


3.5  *  Wavelet  Reconstruction  of  Signals  and  Images  cor¬ 

rupted  with  non-Gaussian  Blur 

A  one-dimensional  signal  5(x)  is  blurred  by  convolution  with  a  filter  F(x)  that  is 
not  a  gaussian  (in  other  words,  the  blurring  process  is  not  equivalent  to  a  diffusion 
in  time,  but  rather  the  result  of  a  more  complicated  process).  The  resulting  blurred 
signal  is  denoted  as  I{x): 

5(x)  *  F(x)  =  7(x)  (199) 

The  goal  is  to  reconstruct  S(x)  by  “antiblurring”  /(x),  which  is  equivalent  to  con¬ 
volving  /(x)  with  the  convolution  inverse  of  the  filter  function  F(x).  In  other  words, 
/(x)  will  be  anti-blurred,  and  the  original  signal  S(x)  will  be  recovered  by  inverting 
Eq.(  199).  To  this  end,  we  will  break  down  the  signal,  the  filter,  and  the  blurred  signal 
into  wavelet  “building  blocks”,  and  process  the  building  blocks  themselves,  instead 
of  working  on  the  entire  waveforms.  We  will  then  prove  that  the  resulting  wavelet 
equations  can  be  easily  inverted,  and  that  the  inverse  problem  is  well-posed,  with 
only  a  mild  restriction  on  the  filtering  function  F(x).  We  proceed  by 


•  expanding  the  original  5(x)  in  terms  of  a  truncated  HR  wavelet  expansion  of 
weight  A: 

S(x)  «  £  flnWnA(x)  (200) 


n=0 


where  {an}n  are  the  Hermite  moments  of  the  signal  5(x),  weighted  by  A. 
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expanding  the  filter  F(x)  in  terms  of  a  truncated  HR  wavelet  expansion  of 
weight  //: 


p(*)  ~  e  ww 

n=0 

where  {6n}„  are  the  Hermite  moments  of  the  filter  F( i),  weighted  by  /z. 


(201) 


•  expanding  the  blurred  signal  I(x)  in  terms  of  a  truncated  HR  wavelet  expansion 
of  weight  v: 

I(x)  «  •£  cjv^x)  (202) 

n= 0 

where  {c„}n  are  the  Hermite  moments  of  the  signal  /(x),  weighted  by  v. 


Eq.(  199)  becomes: 


*  {e^m}  =  E^'M  (203) 


From  Appendix  D,  we  know  that 


Em  l z/^a2+,‘J  ( 


W*(x)  *  W-(x)  =  (x) 


where 


— 

/n,m 


f  (n  +  m)! 


(204) 


(205) 


•  n,m  ,  .  0  njrn  W  _  I  I  V ~  / 

(A2+//2)  2-  V  n!m! 

and  thus,  writing  the  subscript  n  m  as  nm,  and  suppressing  the  superscript  A,(1  when 
confusion  does  not  arise,  Eq.(  203)  becomes: 

Co  =  7oogcA) 

Ci  =  7oiao^i  +  7ioGi^o 

c2  =  702a0^2  +  7llal^l  +  720^2^0 


Cn  —  ^  "  7p.n— pGp^n— f 

P= 0 


(206) 


In  matrix  notation, 

^  Co  \  /  7oo6o 

Ci  7oi^i 

c2  702&2 


0 

7io^o 


0 

0 

720^0 


(207) 


TOn^n  Tfl.n  —  1  ^n  —  1  *Y2,n  — 2^n  — 2  *  •  •  T^nO^O  /  V 
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It  can  be  seen  that  the  above  matrix  is  lower  triangular,  and  thus  that  the  determinant 
is  equal  to  the  product  of  the  main  diagonal,  which  is  equal  to  b „  n*=o  7*o-  That 
pioduct  is  always  nonzero,  provided  that  bo  be  non-zero.  In  other  words,  the  matrix 
equation  (  207)  can  always  be  inverted,  provided  that  the  zeroth  moment  bo  of  the 
filter  be  non-zero.  But  since  we  require  that  the  filter  conserve  the  total  intensity  of 
the  original  waveform  5(x),  bo  must  be  equal  to  1  (see  Appendix  F),  and  Eq.(  207) 
can  always  be  inverted.  The  solution  is  obtained  by  simple  Gaussian  elimination  with 
the  following  recursion  formula: 

1 

Uo  =  - j— Co 

7ooOo 

l  (  »-i  1 

Q n  —  7~  1  ^  ^  7p.n-pap^n-p  f  (208) 

7n0  bo  (  p_0  J 

Once  the  coefficients  {a„}„  are  determined,  the  original  signal  S(x)  can  be  recon¬ 
structed  by  its  wavelet  expansion: 

S(r)  =  !>„»'*(*)  (209) 

n 

Experiments  involving  the  reconstruction  of  simulated  landscapes  and  nebulae 
corrupted  with  non-gaussian  blur  are  being  developed. 
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(a)  IMAGE  FORMING  SYSTEM 


scene  radiance 
f(x.y) 


optical 

system 


2 


image  irradiance 
C(x.y) 


(b)  IMPULSE  RESPONSE 


point  spread 
function 


(c)  LINEARITY 


(d)  SHIFT-INVARIANCE 


1.  Properties  of  linear  and  time-invariant  (LTI)  systems. 
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with  a  smooth  point- ^unction  is  a  “smoothing”  operation. 
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4.  Approximations  of  the  convolution  inverse  of  a  gaussian  filter  of  increasing  width 
f,  as  given  by  Fourier  methods.  The  limits  of  the  Fourier  integral  are  cut-off 
at  ±a:  F'(a,z,|)  =  ^  /*“  e"'*x+V  dk.  The  skirt  in  each  figure  ends  at  the 
minimum  of  F(a.o,x,  |),  while  the  upper  endpoint  of  the  vertical  F  axis  ends  at 
its  maximum. 


46 


HR  Wavelet,  order  0  HR  Wavelet,  order  1  HR  Wavelet,  order  2  HR  Wavelet,  order  3 


HR  Wavelet,  order  4  HR  Wavelet,  order  5  HR  Wavelet,  order  6  HR  Wavelet,  order 


HR  Wavelet,  order  8  ’  HR  Wavelet,  order  9  HR  Wavelet,  order  10  HR  Wavelet,  order  11 


5.  Elements  of  the  sequence  of  ID  HR  wavelets  of  weight  A  =  2,  evaluated  in  the 
range  [-6,6],  The  HR  wavelets  are  rescaled  to  unit  height. 
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HR  W  avelet:  order  12,12  HR  W  avelet:  order  13,13  H  R  Wavelet:  order  1 4,1 4  HR  W  avelet:  order  1  5.1 : 


6.  Elements  of  the  sequence  of  2D  HR  wavelets  of  weight  A  =  /z  =  2,  evaluated 
in  the  range  [-6,6] x [-6,6],  with  the  ordinate  in  the  range  [0,1]  exhibited  in  grey 
scales  (i.e.,  negative  intensities  are  not  plotted).  The  HR  wavelets  are  rescaled 
to  unit  maximum  intensity. 
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7.  First  16  elements  of  a  family  of  80  artificial  “landscapes”.  The  landscapes  were 
obtained  by  summing  five  gaussians  whose  widths  and  offsets  were  obtained 
by  a  random-number  generator  in  an  appropriate  range.  The  landscapes  are 
subsequently  rescaled  to  unit  height. 


ID  HR  wavelet,  order  10000 


8.  ID  HR  wavelets  of  order  1000  and  10000  respectively,  and  of  weight  A  =  1.  The 
renormalization  of  Hermite  polynomials  to  Hd  polynomials  allows  the  numer¬ 
ical  evaluation  of  HR  wavelets  of  high  order  by  eliminating  the  need  for  any 
evaluation  of  factorial  products. 


50 


10.  The  “futuristic”  landscape  68  is  approximated  by  an  HR  wavelet  expansion  of 
width  A  =  1.  The  partial  sums  are  plotted  in  increments  of  20.  The  capturing  of 
very  localized  features  (i.e  high  frequency  content)  is  slower  than  the  capturing 
of  global  behavior  (i.e.  low  frequency  content  features). 
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11.  The  “futuristic”-looking  landscape  71  is  approximated  by  an  HR  wavelet  ex¬ 
pansion  of  width  A  =  1.  The  partial  sums  are  plotted  in  increments  of  20.  The 
wavelet  expansion  experiences  difficulty  in  approximating  one  specific  feature 
of  the  landscape.  This  figure  is  to  be  compared  with  the  next  figure. 


53 


12.  The  “futuristic’Mooking  landscape  71  is  approximated  by  an  HR  wavelet  ex¬ 
pansion  of  width  A  =  1.15.  The  partial  sums  are  plotted  in  increments  of  20. 
This  wavelet  expansion,  as  opposed  to  the  previous  expansion  for  which  A  =  1, 
now  quickly  captures  all  features  of  the  landscape.  The  text  explains  why  this 
happens. 
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14.  The  “realistic’Mooking  landscape  11  is  approximated  by  an  HR  wavelet  expan¬ 
sion  of  width  A  =  1.0.  The  partial  sums  are  plotted  in  increments  of  1.  The  first 
9  partial  sums  are  entirely  undifferentiated,  and  look  exactly  like  the  sequence 
of  underlying  HR  wavelets.  The  approximation  starts  taking  shape  with  the 
13<H  partial  sum. 
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15.  The  “realistic”-looking  landscape  11  is  approximated  by  an  HR  wavelet  expan¬ 
sion  of  width  A  =  1.48.  The  partial  sums  are  plotted  in  increments  of  1.  Ihe 
features  of  landscape  11  are  generated  with  very  few  partial  sums.  It  takes  only 
7  partial  sums  to  recognize  the  landscape.  The  weight  A  =  1.48  of  the  wavelet 
expansion  is  optimal  for  this  landscape. 


16.  The  “realistic’Mooking  landscape  11  is  approximated  by  an  HR  wavelet  expan¬ 
sion  of  width  A  =  2.0.  The  partial  sums  are  plotted  in  increments  of  1.  A  is 
once  again  in  a  suboptimal  range,  as  it  takes  at  least  twelve  partial  sums  to 
recognize  the  landscape. 
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18.  The  HR  moments  for  landscape  11  were  quantized  to  a  range  of  256  =  2R  values. 
In  other  words,  the  values  of  the  HR  moments  were  slightly  modified  to  fall  into 
256  equally  spaced  bins  in  the  range  of  the  HR  moments.  The  weight  A  of  the 
expansion  is  1.48  (optimal).  The  rate  of  convergence  of  the  quantized  wavelet 
expansion  is  only  slightly  modified.  The  HR  wavelet  expansion  is  robust  with 
respect  to  small  errors  on  the  HR  moments. 

60 


19.  First  16  elements  of  a  family  of  80  artificial  “nebulae”.  The  nebulae  were  ob¬ 
tained  by  summing  five  2D  gaussians  whose  widths  and  offsets  were  obtained 
by  a  random-number  generator  in  an  appropriate  range.  The  nebulae  are  sub¬ 
sequently  rescaled  to  unit  maximum  intensity. 
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20.  Letter  “T”  (for  Tina),  which  has  been  diffused  for  a  duration  of  time  t  — 
24,  corresponding  to  12x10*  numerical  iterations  of  the  discretized  diffusion 
equation  (with  unit  diffusivity),  and  subsequently  reconstructed  back  to  t  =  0 
(with  intermediary  results)  by  the  same  number  of  iterations  of  the  discretized 
antidiffusion  equation  (with  unit  diffusivity).  The  reconstruction  is  successful. 
The  text  gives  the  correct  scale  of  time  t. 


letter  T,  t=0 


diffused.  t=30 


i - r  r 


restored,  t=25  restored,  t=20 


restored,  t  =  15 


restored,  t=5 


restored,  t  =  10 


restored,  t  =  0 


21.  Letter  “T”,  which  has  been  diffused  for  a  duration  of  time  t  =  30  (15xl04  itera¬ 
tions  of  the  discretized  diffusion  equation),  and  subsequently  reconstructed  back 
to  t  =  0  (with  intermediary  results)  by  iterating  the  discretized  antidiffusion 
equation  (15xl04  iterations).  The  t  =  0  reconstruction  is  not  successful. 


IIP 

::Sj 
III! 


'  I 


letter  E,  t=0 


diffused,  t=25.2 


:::::: 

...... 


j'llli! 


restored,  t=21.0  restored,  t=16.8 


pi 

mm  1 

lily  :: 

pi 

I  I 

§| 

JSilS 

restored,  t=12.6  restored,  t=8.4 


restored,  t=4.2  restored,  t=0.0 


22.  Letter  “E”  (for  Entropy),  which  has  been  diffused  for  a  duration  of  time  t  =  24 
(12x  104  iterations  of  the  discretized  diffusion  equation),  and  subsequently  re¬ 
constructed  back  to  t  =  0  (with  intermediary  results)  by  iterating  the  discretized 
anti  diffusion  equation  (12xl04  iterations).  The  reconstruction  is  successful. 
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letter  E,  t=0  diffused,  t=30 


restored,  t=25  restored,  t=20 


restored,  t=15  restored,  t=10 


restored,  t=5  restored,  t  =  0 


23.  Letter  r,  ,  which  has  been  diffused  for  a  duration  of  time  t  =  30  (15x10'’  itera¬ 
tions  of  the  discretized  diffusion  equation),  and  subsequently  reconstructed  back 
to  t  —  0  (with  intermediary  results)  by  iterating  the  discretized  antidiffusion 
equation  (15xl04  iterations).  The  t  =  0  reconstruction  is  not  successful. 
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spiral  pattern,  t=0  diffused,  t  =  2 


diffused,  t=24  restored,  t= 


diffused,  t=25.8  restored,  t  =  0 


24.  Spiral  pattern,  which  is  diffused  for  durations  t=24,  25.8,  26.4,  and  27,  and  sub¬ 
sequently  reconstructed  back  to  t  —  0  by  iterating  the  discretized  antidiffusion 
equation.  As  time  is  increased,  the  reconstruction  is  increasingly  unsuccessful. 


+  noise  (SNR=102) 


restored,  t=0 


4-  noise  (SNR=10:>)  restored,  t= 


0 


25.  Letter  “T”,  which  is  diffused  for  t= 12,  and  corrupted  by  additive  noise  at  SNR’s 
=  1,  10  ,  10s,  106,  and  107.  The  reconstruction  by  iterating  the  discretized 
antidiffusion  equation  is  subsequently  attempted.  The  reconstruction  is  only 
successful  for  extremely  high  SNR’s. 


original  signal 


|  pension  of  the  dif- 

i  fused  and  noisy  sig- 

|  no!,  with  weight  A  - 

3.  The  reconstruction 
i  is  successful. 


diffused  signal  for 
t  =  0.3 


diffused  signal  for 
t=0.3,  corrupted  with 
3%  additive  noise  I 


reconstructed  signal, 
antidiffusing  the  50''' 
partial  sum  HR  ex¬ 
pansion  of  the  dif-  I 

fused  and  noisy  sig-  | 

nnl,  with  weight  A  =  j 

3.  The  reconstruction  | 

is  not  successful.  ' 


26.  Landscape  11,  which  is  diffused  for  t  =  0.3  by  convolution  with  the  Green’s 
function  of  the  1+1  diffusion  equation.  The  diffused  signal  is  subsequently 
corrupted  with  3%  additive  noise  (the  noise  field  is  added  uniformly  onto  the 
diffused  signal,  at  a  constant  proportion  of  3%  of  the  maximum  intensity  of  the 
signal).  The  reconstruction  by  antidiffusion  of  the  HR  wavelet  expansion  of  the 
diffused  and  noisy  signal  is  attempted.  The  weight  of  the  HR  wavelet  expansion 
is  A  =  3.  The  25<h  partial  sum  of  the  reconstruction  is  a  good  approximation 
to  the  original  signal.  The  50th  partial  sum  is  not.  The  text  explains  why  this 
happens. 
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original  imagr 


diffused  image  for 
t  =  O.G,  enhanced  for 
contrast 


noise  field 


i 

I 


i 


i  t _ _ _ 

reconstructed  image, 
antidiffusing  the  21" 
partial  sum  HR  ex¬ 
pansion  of  the  dif¬ 
fused  and  noisy  im¬ 
age.  with  weight  A  = 
2.  The  reconstruction 
is  successful. 


diffused  image 

for  t  =  0.6,  corrupted 
with  30%  multiplica¬ 
tive  noise,  enhanced 
for  contrast 


reconstructed  image, 
antidiffusing  the  20  ' 
partial  sum  HR  ex¬ 
pansion  of  the  dif¬ 
fused  and  noisy  im¬ 
age,  with  weight  A  = 
2.  The  reconstruction 
is  not  successful. 


27.  Nebula  62,  which  is  diffused  for  t  =  0.6  by  convolution  with  the  Green’s  function 
of  the  2+1  diffusion  equation.  The  diffused  signal  is  subsequently  corrupted 
with  30%  multiplicative  noise  (the  noise  field  is  added  pixel-by-pixel  onto  the 
diffused  image,  at  a  varying  proportion  of  30%  of  the  intensity  of  the  pixel). 
The  reconstruction  by  antidiffusion  of  the  HR  wavelet  expansion  of  the  diffused 
and  noisy  image  is  attempted.  The  weights  of  the  2D  HR  wavelet  expansion 

n+m  =  21 

are  A  =  fi  =  2.  The  21'*  partial  sum  (ie.  £„£„,)  of  the  reconstruction  is  a 
good  approximation  to  the  original  image.  The  29^  partial  sum  is  not.  The 
text  explains  why  this  happens. 
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28.  Positivity  and  total  intensity  conservation  test  for  the  reconstructions  in  the  two 
previous  figures.  These  tests  yield  the  order  of  the  best  possible  reconstruction 
by  wavelet  antidiffusion  of  blurred  waveforms.  As  these  tests  fail  (as  curves 
2  and  3  diverge),  the  reconstruction  fails  as  well  (the  least-squares  fit  of  the 
reconstruction  with  the  original  waveform,  i.e.  curve  1,  diverges). 


70 


diffused  image  for 
t  =  0.2,  enhanced  for 
cont  rast 


diffused  image  for 


t=0.2 


discretization  of  the 
antidiffusion  equation 
with  ' V H  =  |04. 


■gl 

Reconstruction  by- 

discretization  of  the 
antidiffusion  equation 
with  5 A  11  -  io,;. 

Reconstruction  by 

discretization  of  the 
antidiffusion  equation 
with  .C.Y H  -  j  x  I0'=. 

-  * 

m 

Reconstruction  by 

discretization  of  the 
antidiffusion  equation 
with  >\V/i r  ~  J0H. 

Reconstruction  bj 

discretization  of  the 
antidiffusion  equation 

i'.V/i  =  10,s  (equiva- 

lent  to  double  preci¬ 
sion  round-off  error) 


29.  Attempt  for  reconstruction  of  the  blurred  image  of  nebula  62,  by  iterating  the 
discretized  antidiffusion  equation.  While  HR  wavelet  antidiffusion  allows  recon 
struction  with  very  ^ow  SNR’s,  it  can  be  seen  that  antidiffusion  by  discretization 
of  the  antidiffusion  equation  is  only  possible  for  very  high  SNR’s. 
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Appendix  A:  Hyperdistribution  Expansion  of  the 
Convolution  Inverse  of  a  Gaussian  Filter 


The  Gaussian  filter  6\(x,y)  is  a  “radial”  filter,  that  is,  it  is  separable  in  its  two 
independent  variables  x  and  y.  As  a  result,  its  hyperdistribution  expansion  is  equally 
separable,  and  the  2D  Bochner  Algebra  reduces  to  the  product  of  two  ID  Algebras  (see 
Eq.(165)).  The  hyperdistribution  expansion  of  its  inverse  is  thus  equally  separable, 
and  We  can  restrict  the  study  to  just  one  variable,  without  loss  of  generality.  The 
hyperdistribution  expansion  in  x  is  obtained  from  Eq.(163),  which  is  repeated  here 

+  oo  \  4  1 

m*)  =  (210) 

„=o  o  n! 

The  ID  Bochner  algebra  thus  reduces  to  determining  the  coefficients  bk  such  that 

n  /  A4  \n  —  k 

{m) 


where  <5n0  =  1  if  n  —  0  and  0  otherwise.  Trying  the  substitution 


in  the  sum  in  Eq.(211)  yields: 


bk  = 


(-T)* 

k\ 


\4  n 

(j>"£ 

°  k=o 


(~l)k 
( n  -  ky.kl 


The  partial  sum  above  can  be  rewritten  as 


f  orw 

to  {n-k)\k\ 


(212) 


(213) 


(214) 


and  the  binomial  expansion  of  (1  —  l)"/n!  is  recognized,  which  is  identically  equal  to 
0,  except  for  when  n  =  0,  in  which  case  it  is  equal  to  1.  q.c.d. 
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Appendix  B:  Discretization  of  the  operator  V2 


In  this  appendix,  we  discretize  the  continuous  operator  V2,  in  order  to  discretize  both 
diffusion  and  antidiffusion  equations.  The  discretization  scheme  will  be  centrat-in- 
space,  and  forward  in  time.  Let  xoo, Xoi^ Xo„, Xio, Xii,...,xi„, ...  denote  the  spatial 
location  of  the  grid  points,  and  let  the  real  function  /,  of  the  variable  xpq,  be  the 
measure  of  a  real- valued  distribution  assigned  to  the  grid.  If  confusion  does  not  arise, 
we  will  denote  f{xpg)  as  fM.  A  Taylor  series  of  order  2x2,  for  ftJ ,  around  xpq  can 
be  written  as: 


/,;  =  EE—  (*•«  -  X^)n  ^  -  X”^m  7nr/nm  + 


(215) 


n=0  m— 0 


If  {xij}ij  are  the  locations  of  all  eight  neighbors  surrounding  xpq,  and  a,:  is  the  weight 
assigned  to  the  neighbor  at  x,j,  we  write: 


3  3 

XIXIQo7u 

.=i j=i 


3  3  2  2 


—  XIXIa,-7  XI  XI  I  1  (I'«  Xnq)  VXj9  Xm,)  fnm  +0(6;) 


=  \  n=0  m=0 


=  XI  XU  ( XT  XI  (*•»  -  xn,)n  (xJq  -  xmqy 

n=0  m=0  ^  i— 0  j— 0 

+  0(61) 


n\m\  dndv 


(216) 


We  are  interested  in  evaluating  V2/p,  as  a  linear  combination  of  the  type  ^  Qi jf,}, 

where  {x.j},^  are  the  locations  of  ali  eight  neighbors  of  xpq.  Since 

d2  d2 

V2/o  =  ^/o  +  ^/o  (21 -) 

We  set  {  }n,m  =  {  }m,„  =  0  for  (m,n)  in  {(0,0),  (0,1),  (1,1)},  and  {  }„,m  = 

{  }m.n  for  (m,n)  =  {(0,2)},  in  Eq.(2I6).  These  constraints  generate  five  equations, 
which  reduce  to  the  following  two  independent  equations: 


Q  £  f  lo  +0  Hf\a  +7/p. 

(a  +  ;m 

-4(a  -f  /?) 


+  wi) 


(218) 

(219) 


where  a  with  no  subscripts  denotes  the  coefficient  assigned  to  each  cell  of  the  diamond 
stencil,  while  0  denotes  the  coefficient  assigned  to  each  cell  of  the  square  stencil,  and 
7  is  the  coefficient  assigned  to  the  center-cell  xpq.  a  and  8  remain  free-  parameters  in 
file  discretization. 


Appendix  C:  Global  Invariance  of  the  Set  of  Para¬ 
metric  Hermite-Rodriguez  Wavelets  in  the  Diffu¬ 
sion  Group 


In  this  appendix,  we  show  that  the  diffusion  or  antidiffusion  of  an  HR  wavelet  yields  a 
rescaled  HR  wavelet  of  same  order,  and  of  respectively  bigger  or  lesser  weight.  Since 
the  Dirac  delta  is  the  unit  element  for  the  operation  of  convolution,  and  since  the 
diffusion  of  the  Dirac  delta  yields  a  Gaussian  of  width  two  times  the  square  root  of 
the  diffusion  duration,  we  can  write  that: 


ll'„x(x)  =  e,v2«x)  .  W‘(*))  =  (e‘vi«(x))  *  VVA(x)  = 


e-W« 
2  s/nt 


*  Wnx(x)  (220) 


The  last  operation  of  convolution  above  can  be  evaluated  explicitly  by  employing  the 
Rodriguez  formula  for  HR  wavelets,  and  the  chain  rule  on  V”: 


e-x3/4< 

2\/tU 


•  tV„A(x)  =  / 


+oo  e-V2/4t 
-oo  2\fxi 


(— 1  )"AnV”_y 


v/tFA  dy 


(— l)r’A”Vn 

f+°°  e~yJ/4t  e~ 

-(x-v)2 

— oo 

y/n\ 

(  — l)”AnVn 

j2 

e 

t/2  (n!) 

\A(^2  +  4t) 

An 

iyVTO7(x) 

(\/A2  -f  4<)T1 

dy 


(221) 

(222) 

(223) 

(224) 


And  consequently,  we  can  write: 

diffusion:  e+,v’lVnA‘(x)  =  ({j)V*»(i)  (225) 

antidiffusion:  e~tV2  W*2(x)  =  [j*}  Hy^x)  (226) 

with  A2  -  \jx\  -f  At  (227) 

In  other  words,  the  set  of  parametric  HR  wavelets  is  a  globally  invariant  subspace 
-in  the  weight  parameter  A-  of  the  diffusion/antidiffusion  group[14].  This  property 
of  HR  wavelets  enables  one  to  formally  diffuse/antidiffuse  one-dimensional  signals, 
or  two-dimensional  images,  by  diffusing/antidiffusing  the  HR  wavelets,  which  are  the 
building  blocks  of  the  signals’  or  images’  HR  expansions. 


*  Appendix  D:  Global  Invariance  of  the  Set  of 
Parametric  Hermite-Rodriguez  Wavelets  in  the 
Convolution  Group 


This  appendix  is  a  generalization  of  the  previous  appendix:  We  show  that  the  con¬ 
volution  of  two  HR  wavelets  of  respective  orders  n  and  m,  and  respective  weights  A 
and  /*,  yields  another  HR  wavelet,  of  order  n  -f  m,  and  weight  y/\2  +  ft2. 

Lemma  1.  The  Rodriguez  formula  for  Hermite  polynomials  yields: 

(-l)"//„(x)e-*2  =  V"(e-*2)  (228) 


It  can  be  rescaled  to  yield: 


(-l)nHn(x/\)e~x'/x2  =  Vn(e-r2/Ai)  □ 


(229) 


Lemma  2.  By  definition  of  HR  v/avelets, 


Wx(x)  =  Hdx(x) 


Hn(x/\)  e-*2'* 
2  uI2\/ti\ 


By  Lemma  1, 


f'_nr>  A*  e-x2l*2 

Wx(x)  =  - _ _  Vn( - ) 

n[  ’  2»/2v/^i  1  v^FA  ’ 


and  thus 


f_1\n+m\n  m  -x‘  / 

H^(x)  .  W£  (X)  =  V"(-=-)  .  V"(~-=-)  □ 

2 2  Vnlml  V^A 


(230) 


(231) 


(232) 


Lemma  3.  The  convolution  of  the  n</l  derivative  of  e  z2^2  with  the  m</l  derivative  of 
e'x  /(1  is  given  by 

Vn(e-x>/A2)  *  Vm(e-l2/"J)  =  f+<X>V^{e-y2,x2)\'^_v{e-^-v)2^2)dy  (233) 

J —oo 

By  the  chain  rule  on  V”_y, 

Vn(e-r2/A2)  *  Vm(e~x2^2)  =  [+°°V”(e-y2/x2)V™(e-{l-y)2/"2)dy 

J  —  OO 

=  vm  r°°  Vny{e-y2/x2)e-(l-y]2lil2dy  (231) 

J  —  OO 

By  a  change  of  variables,  and  the  chain  rule  on  V"_y, 

Vn(c'rVA2)  *  Vm(c-x2/‘l2)  =  Vn+m  /+°°  c~y2/x2  c-[r-y)2/,,2dy  (235) 

J  —  OO 
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and  finally, 


yn(e-*2/AJ)  +  = 


A2  -f  (i‘ 


yn+m(e-i1/(iJ+|i!))  □ 


(236) 


Finally  combining  Lemmas  2  and  3, 


Wfa)  *  W£(Z)  =  7n,m  (*) 


(237) 


In  other  words,  the  set  of  parametric  HR  wavelets  is  a  globally  invariant  subspace  -in 
the  weight  parameter  A-  under  convolution.  This  property  of  HR  wavelets  enables  one 
to  formally  convolve/deconvolve  one-dimensional  signals,  or  two-dimensional  images, 
by  convolving/deconvolving  the  HR  wavelets,  which  are  the  building  blocks  of  the 
signals’  or  images’  HR  expansions. 
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Appendix  E:  Waveform  Total  Intensity  Conserva¬ 
tion  in  the  Diffusion  Group 

We  first  prove  that  any  diffusion/antidiffusion  process  conserves  the  total  intensity  of 
the  initial  condition.  In  other  words,  if  We  denote  by  I(x)  the  initial  condition  of  a 
diffusion/antidiffusion  process,  then  the  total  intensity  is  given  by  I(x)dx,  and 

f+°°  e±tV'l(x)dx  =  0  (239) 

ut  J — oo 


The  proof  is  as  follows.  By  integrating  in  space  the  diffusion/antidiffusion  equation 
for  u(x,  t): 


r+ oo  r+oo 

j  — -dx±  j  V2udx  =  0 

J  —  oo  Ot  J — oo 

(240) 

With  the  help  of  Gauss’  theorem, 

3  r+°°  r  1+00 

—  /  udxifc  n-Vu  =0 

3t  J- oo  1  J  -oo 

(241) 

and  thus 

d  r+°°  , 

—  /  udx  -  0 

ut  J — oo 

(242) 

Since  u(x,t)  =  e±tV 2  /(x),  we  obtain  the  desired  result: 

Jt  =  0 

(243) 

The  same  reasoning  can  be  easily  extended  to  the  case  where  the  initial  condition  is 
a  two-dimensional  image  T(x,y). 

Furthermore,  the  total  intensity  of  a  truncated  HR  wavelet  expansion  of  any  wave¬ 
form  is  exactly  equal  to  the  total  intensity  of  the  waveform.  That  is  because  the  total 
intensities  of  all  HR  wavelets  are  identically  equal  to  zero,  except  for  the  total  inten¬ 
sity  of  the  HR  wavelet  of  order  0,  which  is  by  definition  the  zeroth  HR  moment,  and 
thus  the  total  intensity  of  the  waveform.  In  other  words,  the  expansion  of  a  waveform 
in  terms  of  a  truncated  HR  wavelet  series  conserves  the  total  intensity  of  the  original 
waveform  (this  property  is  also  true  in  Fourier  series). 


77 


"Appendix  F:  Do  Point-Spread-Functions  Conserve 
the  Total  Intensity  of  the  Original  Waveform? 

In  this  appendix,  we  derive  the  constraint  imposed  on  the  point-spread  function  of 
LTI  systems  which  conserves  the  total  intensity  X  of  their  input.  In  other  words, 
if  u(i)  denotes  the  input  function,  r(i)  denotes  the  output,  and  f(x)  denoted  the 
impulse  response  or  point-spread-function,  we  require  that 


/too  r+oo 

v(x)dx  =  /  u(x)dx  =  X 

-OO  J  —  OO 


(244) 


/+oo  r+oo 

f(x  —  x')u(x')dx'  =  /  f(x')u(x  —  x')dx'  (245) 

-OO  J  —  OO 


and  thus,  by  exchanging  the  order  of  integration, 


r+oo  rr+°°  r+oo  f  r+oo  1 

I  v{x)dx  =  I  f(x')u(x  —  x')dx'dx  =  I  f(x')  |  /  u(x  —  x')dx\  dx' 

J  —  OO  J  J  —  oo  J —OO  W—  OO  J 

/  +  OO 

f{x')Xdx'  (246) 

-oo 


We  conclude  that  the  sufficient  condition  for  total  intensity  conservation  of  the  input 
is  that  the  total  intensity  of  the  point-spread-function  be  normalized  to 


/+oo 

f(x)dx  =  1 

-OO 


(247) 
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